{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "name": "01_Consistent.ipynb", "provenance": [], "toc_visible": true, "include_colab_link": true }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "fWMuDsTQjOCi" }, "source": [ "# Consistency of a Multistep method\n", "\n", "#### John S Butler \n", "john.s.butler@tudublin.ie \n", "[Course Notes](https://johnsbutler.netlify.com/files/Teaching/Numerical_Analysis_for_Differential_Equations.pdf) [Github](https://github.com/john-s-butler-dit/Numerical-Analysis-Python)\n", "\n", "## Overview\n", "A one-step or multistep method is used to approximate the solution of an initial value problem of the form\n", "\\begin{equation} \\frac{dy}{dt}=f(t,y), \\end{equation} \n", "with the initial condition\n", "\\begin{equation} y(a)=\\alpha.\\end{equation} \n", "The method should only be used if it satisfies the three criteria:\n", "1. that difference equation is __consistent__ with the differential equation;\n", "2. that the numerical solution is __convergent__ to the exact answer of the differential equation;\n", "3. that the numerical solution is __stable__.\n", "\n", "In the notebooks in this folder we will illustate examples of consisten and inconsistent, convergent and non-convergent, and stable and unstable methods. \n", "\n", "## Introduction to Consistency\n", "In this notebook we will illustate an __inconsistent__ method." ] }, { "cell_type": "code", "metadata": { "id": "ObAwcaFNQpOf", "outputId": "d73494fb-88ce-4b6b-f1de-88ccf90a763e" }, "source": [ "from IPython.display import HTML\n", "HTML('')" ], "execution_count": null, "outputs": [ { "output_type": "stream", "text": [ "/Users/johnbutler/anaconda3/lib/python3.7/site-packages/IPython/core/display.py:717: UserWarning: Consider using IPython.display.IFrame instead\n", " warnings.warn(\"Consider using IPython.display.IFrame instead\")\n" ], "name": "stderr" }, { "output_type": "execute_result", "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": { "tags": [] }, "execution_count": 1 } ] }, { "cell_type": "markdown", "metadata": { "id": "mIJEJ0OUQpOm" }, "source": [ "### Definition\n", "A one-step and multi-step methods with local truncation error $\\tau_{i}(h)$ at the $i$th step is said\n", "to be __consistent__ with the differential equation it approximates if \n", "\\begin{equation}\\lim_{h \\rightarrow 0} (\\max_{1 \\leq i \\leq N}|\\tau_{i}(h)|)=0 \\end{equation} \n", "where\n", "\\begin{equation}\\tau_{i}(h)=\\frac{y_{i+1}-y_{i}}{h}-F(t_i,y_i,h,f) \\end{equation} \n", "As $h \\rightarrow 0$ does $F(t_i,y_i,h,f) \\rightarrow f(t,y)$. \n", "\n", "All the Runge Kutta, and Adams methods are consistent in this course. This notebook will illustrate a non-consistent method which with great hubris I will call the Abysmal-Butler methods.\n", "\n", "## 2-step Abysmal Butler Method \n", "\n", "The 2-step Abysmal Butler difference equation is given by\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{2}(4f(t_i,w_i)-3f(t_{i-1},w_{i-1})), \\end{equation} \n", "which can be written as \n", "\\begin{equation}\\frac{w_{i+1} -w_{i}}{h} = \\frac{1}{2}(4f(t_i,w_i)-3f(t_{i-1},w_{i-1})). \\end{equation} \n", "\n", "## Intial Value Problem\n", "To illustrate consistency we will apply the method to a linear intial value problem given by\n", "\\begin{equation} y^{'}=t-y, \\ \\ (0 \\leq t \\leq 2),\\end{equation} \n", "with the initial condition\n", "\\begin{equation}y(0)=1,\\end{equation} \n", "with the exact solution\n", "\\begin{equation}y(t)= 2e^{-t}+t-1.\\end{equation} " ] }, { "cell_type": "markdown", "metadata": { "id": "sOosaQngjOCj" }, "source": [ "## Python Libraries" ] }, { "cell_type": "code", "metadata": { "id": "TXMyUzxojOCk" }, "source": [ "import numpy as np\n", "import math \n", "import pandas as pd\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "GovXdrE0jOCo" }, "source": [ "### Defining the function\n", "$$ f(t,y)=t-y.\\end{equation} " ] }, { "cell_type": "code", "metadata": { "id": "pAx9mdLUjOCp" }, "source": [ "def myfun_ty(t,y):\n", " return t-y" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "OhMnF7tgjOCs" }, "source": [ "## Discrete Interval\n", "Defining the step size $h$ from the interval range $a \\leq t\\leq b$ and number of steps $N$ \n", "\\begin{equation}h=\\frac{b - a}{N}.\\end{equation} \n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=t_0+ih,\\end{equation} \n", "where $t_0=a.$\n", "\n", "Here the interval is $0\\leq t \\leq 2$ and number of step 40 \n", "\\begin{equation}h=\\frac{2−0}{40}=0.05.\\end{equation} \n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=0+i0.5,\\end{equation}\n", "for $i=0,1,⋯,40.$" ] }, { "cell_type": "code", "metadata": { "id": "LOnnFiILjOCt", "outputId": "ab5227bd-6d3b-49c6-b100-53f53651d228" }, "source": [ "# Start and end of interval\n", "b=2\n", "a=0\n", "# Step size\n", "N=40\n", "h=(b-a)/(N)\n", "t=np.arange(a,b+h,h)\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(t,0*t,'o:',color='red')\n", "plt.xlim((0,2))\n", "plt.title('Illustration of discrete time points for h=%s'%(h))\n", "plt.show()" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAEICAYAAAAN2sPKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAee0lEQVR4nO3de7glVX3m8e8rCA4XucjFbrqhUUkiJmqgQzBmREN0AEXUJxlxGLxlJMSQkBkTgxpz0RhNmFFCJBJiMCpE4gwaOgYCRGOcQUUb5SJp0IYo3dJIc7+piP7mj6qT3r1773N29znSvc75fp5nP72raq1aq1ZV135P1a5zUlVIkiRp2/aYrd0BSZIkzczQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5sWtCSvTvL/BqYryVO2Zp/GSXJ2krduhXZ/Jcm3kjyQ5AkTlP96kp/v3785yft/+L384UlyQpLLtnY/ppNk/37/bPcotPXSJGv69n5yDtb36ST/bS76Js13hjYtOIOh4oe0/r9O8oezXMdGYRKgqk6uqrfPrneb3Y/HAu8GXlBVu1TVnZtTv6r+qKoetQ/kJMv64L39XNWvqvOr6gVz18u5V1W39Pvn+zOVne0YAf8TOKVv78tbuI45l+SZSa5K8lD/7zOnKbtjknOT3JfktiT/Y2h5JXmwD6YPtP6Dh+YPQ5v0KJvFh+XWsC/wOOD6rd0RaG7s5qsD2MLj4Yd1JTDJDsBFwHnAHsAHgYv6+aP8PnAQ3bY8D3hjkqOGyjyjD6a7PJo/eEjTMbRJYwzfthm8+pXOe5LcnuTeJNcm+fEkJwEn0H0IPJDk7/vyX0/y20muBR5Msn2S05LclOT+JP+a5KV92acCZwPP6tdxTz9/oyt4SV6XZHWSu5KsSLJ4YFklOTnJ15LcneSsJBmznTsmOSPJrf3rjH7ejwA39sXuSfKpMfVPTPKNJHcmecvQst9Pcl7//nFJzuvL3ZPki0n27ZftmeQDfft3J/m7fv5zk6ztx+424ANJHjMwdncm+WiSPfsmPzPQ3weSPKtfz2uTrOrXfWmSA8bs9k3qD1/17Mf29f3Y3p/k7UmenORz/ZWbjw6GhSQvSnJ1v82fTfL0MW1PrfvXk9yc5I4kpyd5TL/sMUl+px/r25N8KMlu/bKNrp71x+7bk1zR9/GyJHtNs41PSfIv/bF8R5K/HdG3HZM8AGwHXJPkpn7+U/v27klyfZIXD9T56yTvS3JxkgfpAtIoB4zp66SeC2wPnFFV362qM4EAPzem/CuBt1fV3VW1CvhL4NWb2ab0qDO0SVvmBcBzgB8BdgdeDtxZVecA5wN/0v+EfuxAnVcALwR2r6pHgJuA/wjsBvwBcF6SRf2HyMnA5/p17D7ceJKfA94J/GdgEfAN4IKhYi8Cfgp4Rl/uP43ZlrcAhwPP7MseBvxOVX0VeFpfZveq2uQDMMnBwPuAE4HFwBOAJWPaeVW/rUv7cicD3+6XfRjYqW9vH+A9A/WeCOxJd1XkJODXgZcAR/Rt3g2c1Zd9zkB/d6mqzyV5CfBm4GXA3sD/BT4ypo+b1B9T7ijgULpxeyNwDl1YXwr8ON2+JskhwLnAL/fb/BfAiiQ7jlkvwEuB5cAhwHHAa/v5r+5fzwOeBOwCvHea9fwX4DV047kD8JvTbOPbgcvorlItAf5seGV9GNqln3xGVT053e3zv+/r7gP8GnB+kh8d6sc7gF2BjW75T9BX+jA47nVaX+xpwLW18d9lvJYNxy8D69uD7ri5ZmD2NSPKfibdrdOPJVk2pt/So8rQJm2Z79F9CP0YkKpaVVXrZqhzZlWtqapvA1TV/66qW6vqB1X1t8DX6ALTJE4Azq2qL1XVd4E30V2ZWzZQ5l1VdU9V3QL8M10oG7eut1XV7VW1ni5AnjhhP34B+ERVfabvx1uBH4wp+z264PKUqvp+VV1VVfclWQQcDZzcX/n4XlX9y0C9HwC/14eGb9MFoLdU1dq+zd8HfiHjb53+MvDOfh89AvwR8MxprrZN4o+r6r6quh74CnBZVd1cVfcClwBTX9B/HfAXVXVlv80fBL5LF/amW/dd/X47gz4A0u2nd/ftPEC3z4+fZrs/UFVf7cfso4zf/9DtmwOAxVX1naoaF66GHU4XHt9VVQ9X1aeATwz0GeCiqrqiP86/s7l9rardp3m9qy+2C3Dv0Drvpfs/OmyXgeXjyh4BLKP7/30r8Ilpxll61BjapC3Qfzi9l+4Kz7eSnJPk8TNUWzM4keSVA7fN7qG7QjPpbaHFdFfXpvrzAHAnsN9AmdsG3j/Ehg+radfVv188puyouv++XVX1YN+PUT4MXApc0N8G/ZP+Ss1S4K6quntMvfVDH/YHAB8fGLdVwPfpvn83ygHAnw6Uv4vu1tl+Y8pP4lsD7789YnpqrA8A3jB4dYhue6cb38HjZHBfjNpP2zN+uyfd/9BdLQzwhf4W52unKTtoMbCmqgaD+jfYeGzXMLPN6esoDwDD//8eD9w/puzU8pFl+x9CHq6qe4BTgQOBp25mn6Q5Z2iTxnuQ7pbdlCcOLqyqM6vqULrbKj8C/NbUojHr+/f5/VWevwROAZ7Q3wL9Ct0H53TrmHIrXSCYWt/OdFexvjlDvRnXBezfz5vEOroQMtWPnfp+bKK/gvYHVXUw8DN0t29fSfehvmeS3ce0MTwWa4Cjh664PK6qvjmi7FT5Xx4q/x+q6rMTtDVba4B3DLW9U1WNuz0LA+PJxvti1H56hI0D4yQ22caquq2qXldVi+muTP55JvvVN7cCS6e+dzfQr8HjcFZjmg1PcI56vbkvdj3w9GSj720+nREPTPQ/HKyj+yrAlGeMKju0DSO/Eyo9mgxt0nhXAy9LslP/AfZLUwuS/FSSn+6vFD0IfIfuag90H6JPmmHdO9N9EKzv1/cauittU74FLMn4p9/+BnhNul9zsCPdLb8rq+rrm7F9Uz4C/E6SvfsvgP8u3VN4k/g/wIuS/Gzf17cx5ryS5HlJfiLdE4T30d2S+35/W/kSuqCwR5LHJnnOqHX0zgbeMXV7s+/3cf2y9XS3U580VP5NSZ7Wl98tyS+OWfeo+rPxl8DJ/bGSJDsneWGSUbftpvxWPw5L6a7yTD0U8BHgvyc5MMkudPv8b/tbvptjk21M8otJpr6LeDfdsTnjrw8BrqQ7/t/Y77fnAsey6fcrt9jAE5yjXn/UF/t0399f7x+YOKWfP/LhGeBDdMf8Hkl+jO429l8DJHla//9qu36c/xddCF01V9skbSlDmzTee4CH6QLUB+keMJjyeLoP5LvpbgfdSff7qwD+Cji4vx32d6NWXFX/Svdh8Ll+/T8BXDFQ5FN0P/nfluSOEfU/Sff9sQvprho8GTh+SzYS+ENgJd0Xt68DvtTPm1H/na5fpQuR6+jGY+2Y4k+kC3n30X0A/gsbwuGJdCHuBuB24DemafZPgRXAZUnuBz4P/HTfn4fovvR+RT/+h1fVx4E/prstex/dFc2jx2zPJvVnGoPpVNVKukDwXrqxWc3MTyleBFxF90PDP9AdT9A90PBhuqc//43uB4Vf24I+jdrGnwKuTPd06Arg1Kr6twnW9TDwYrrxvAP4c+CVVXXD5vZrNvp+vITuyu09dA9vvKSfP/ULkgevpP0e3YNA36A7Dk+vqn/sl+1LF5TvA26m+27bi6rqez/0DZFmkI0ftpEkbS1JCjioqlZv7b5I2vZ4pU2SJKkBhjZJkqQGeHtUkiSpAV5pkyRJakCTv+F5r732qmXLlm3tbkiSJM3oqquuuqOq9p7tepoMbcuWLWPlypVbuxuSJEkzSvKNmUvNzNujkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDZiT0JbkqCQ3Jlmd5LQRy5PkzH75tUkOGVq+XZIvJ/nEXPRHkiRpvpl1aEuyHXAWcDRwMPCKJAcPFTsaOKh/nQS8b2j5qcCq2fZFkiRpvpqLK22HAaur6uaqehi4ADhuqMxxwIeq83lg9ySLAJIsAV4IvH8O+iJJkjQvzUVo2w9YMzC9tp83aZkzgDcCP5iukSQnJVmZZOX69etn1WFJkqTWzEVoy4h5NUmZJC8Cbq+qq2ZqpKrOqarlVbV877333pJ+SpIkNWsuQttaYOnA9BLg1gnLPBt4cZKv091W/bkk581BnyRJkuaVuQhtXwQOSnJgkh2A44EVQ2VWAK/snyI9HLi3qtZV1ZuqaklVLevrfaqq/usc9EmSJGle2X62K6iqR5KcAlwKbAecW1XXJzm5X342cDFwDLAaeAh4zWzblSRJWkhSNfz1s23f8uXLa+XKlVu7G5IkSTNKclVVLZ/tevyLCJIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1YE5CW5KjktyYZHWS00YsT5Iz++XXJjmkn780yT8nWZXk+iSnzkV/JEmS5ptZh7Yk2wFnAUcDBwOvSHLwULGjgYP610nA+/r5jwBvqKqnAocDvzqiriRJ0oI3F1faDgNWV9XNVfUwcAFw3FCZ44APVefzwO5JFlXVuqr6EkBV3Q+sAvabgz5JkiTNK3MR2vYD1gxMr2XT4DVjmSTLgJ8ErpyDPkmSJM0rcxHaMmJebU6ZJLsAFwK/UVX3jWwkOSnJyiQr169fv8WdlSRJatFchLa1wNKB6SXArZOWSfJYusB2flV9bFwjVXVOVS2vquV77733HHRbkiSpHXMR2r4IHJTkwCQ7AMcDK4bKrABe2T9Fejhwb1WtSxLgr4BVVfXuOeiLJEnSvLT9bFdQVY8kOQW4FNgOOLeqrk9ycr/8bOBi4BhgNfAQ8Jq++rOBE4Hrklzdz3tzVV08235JkiTNJ6ka/vrZtm/58uW1cuXKrd0NSZKkGSW5qqqWz3Y9/kUESZKkBhjaJEmSGmBokyRJaoChTZIkqQGGNkmSpAYY2iRJkhpgaJMkSWqAoU2SJKkBhjZJkqQGGNokSZIaYGiTJElqgKFNkiSpAYY2SZKkBhjaJEmSGmBokyRJaoChTZIkqQGGNkmSpAYY2iRJkhpgaJMkSWqAoU2SJKkBhjZJkqQGGNokSZIaYGiTJElqgKFNkiSpAYY2SZKkBhjaJEmSGmBokyRJaoChTZIkqQGGNkmSpAYY2iRJkhpgaJMkSWrAnIS2JEcluTHJ6iSnjVieJGf2y69NcsikdUe66ipYtgzOP3+yDp5/flf+MY+x3kKp10IfrWc967VXr4U+Wm+bq3coHDpZhRlU1axewHbATcCTgB2Aa4CDh8ocA1wCBDgcuHLSuqNeh0IVVO20U9V559W0zjuvKzdVx3rzv14LfbSe9azXXr0W+mi9bbLeoVA1Q7aZ5DX7FcCzgEsHpt8EvGmozF8ArxiYvhFYNEndUa9DBwdu//2rjjii6sMf7gbpwQe76Qsu6KaXLt14oKdee+/dLV+3rit/ySXd9C23dNP77DO63uLF3fIrrujKX3ddN/2FL3TTixaNrrdoUbf8iiu68jfc0E1/+tPd9OLFo+vts0+3fN26rvyKFd30+vXd9F57ja63dGm3/MMf7so//HA3/YEPdNMHHDC63h57VB111IaD7owzqo49dsP07ruPrnfAAd3yt72t6oQTNpR/61urXv3q8e3tumvV61+/ofypp3avKbvsMn17r35118aUE06o2m238XVe9rKq00/fUP7YY7ttrBrfxz333FD+iCO6MazqxvSII6qe8ITx++CII6ouvLArv359N71iRTe9ZMn4fV5VddNNXflPf7qbvuGGbnrffccfY0ccUfXlL3flv/CFbvq667rpcfUWL+6WX355V/6WW7rpSy7ppsf1c6+9uuX33NOVv+CCbvrBB7vpceMyte/OOafqyCM3jO1ZZ3XH3rj9sPvu3f6b8s53Vr385Rump9vvVVWnnVb1utdtKP+GN3TH3rj2dtmlKzPlda/r1jFl552nb+/lL+/6OGXq2BvX3uMe143BlCOP7MZoyo47jm9v+Lx3zz0bjr1x7e2446bnvcsv76Zvumn69obPe1/+8oZjb7r2hs97N93UTV9++fTtDZ/3Lrxww7E3XXvD570p55zTjfe49obPe6ef3u2/cW3tttvo896Uxz9++mNl+Lz3+td3x9649nbeedPz3tvetmF6OGAMtzd43qvq/t+dddb0Yzl83hv8zJ1u3w2f9wY/c6drb/i8N/iZO117w+e9wc/c6dobPu8NfuZO197weW/wM3egvbkKbXNxe3Q/YM3A9Np+3iRlJqkLQJKTkqxMsnKjBWvWjCo+sMa1o+ffccf09davHz1/3brp69122+bNn2m94/ox5c47R88ft91Tbrll9Py7756+3j33bN76Zlp+//3T13vggc1v7957N7/OdMvvumv6elu6D775zdHzZ9rnt98+ev5Mx9i4ejMd0+P6OW67Z1q+pfth3LE3Za73+7hjb8qDD85te9/5zvT1vvvduW1v3PoWSnvjxnu69sYtG3fsTbnvvs1va7rl4469KQ89NLftbWv7rvX2ZmO2qQ/4ReD9A9MnAn82VOYfgJ8dmP4k3f3dGeuOeh06nHSnMy5ZW2/+1muhj9aznvXaq9dCH623Tdbz9ihs8/ewrbeV6rXQR+tZz3rt1Wuhj9bbJuttS6Fte+Bm4EA2PEzwtKEyL2TjBxG+MGndUa9DpxLuTAM2OHAHHFCVWG+h1Guhj9aznvXaq9dCH623zdWbq9CW6sLTrCQ5BjiD7mnQc6vqHUlO7m+/np0kwHuBo4CHgNdU1cpxdWdqb/ny5bVy5cqZikmSJG11Sa6qquWzXs9chLZHm6FNkiS1Yq5Cm38RQZIkqQGGNkmSpAYY2iRJkhpgaJMkSWqAoU2SJKkBhjZJkqQGGNokSZIaYGiTJElqgKFNkiSpAYY2SZKkBhjaJEmSGmBokyRJaoChTZIkqQGGNkmSpAYY2iRJkhpgaJMkSWqAoU2SJKkBhjZJkqQGGNokSZIaYGiTJElqgKFNkiSpAYY2SZKkBhjaJEmSGmBokyRJaoChTZIkqQGGNkmSpAYY2iRJkhpgaJMkSWqAoU2SJKkBhjZJkqQGzCq0JdkzyeVJvtb/u8eYckcluTHJ6iSnDcw/PckNSa5N8vEku8+mP5IkSfPVbK+0nQZ8sqoOAj7ZT28kyXbAWcDRwMHAK5Ic3C++HPjxqno68FXgTbPsjyRJ0rw029B2HPDB/v0HgZeMKHMYsLqqbq6qh4EL+npU1WVV9Uhf7vPAkln2R5IkaV6abWjbt6rWAfT/7jOizH7AmoHptf28Ya8FLpllfyRJkual7WcqkOSfgCeOWPSWCdvIiHk11MZbgEeA86fpx0nASQD777//hE1LkiTNDzOGtqr6+XHLknwryaKqWpdkEXD7iGJrgaUD00uAWwfW8SrgRcCRVVWMUVXnAOcALF++fGw5SZKk+Wi2t0dXAK/q378KuGhEmS8CByU5MMkOwPF9PZIcBfw28OKqemiWfZEkSZq3Zhva3gU8P8nXgOf30yRZnORigP5Bg1OAS4FVwEer6vq+/nuBXYHLk1yd5OxZ9keSJGlemvH26HSq6k7gyBHzbwWOGZi+GLh4RLmnzKZ9SZKkhcK/iCBJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDZhVaEuyZ5LLk3yt/3ePMeWOSnJjktVJThux/DeTVJK9ZtMfSZKk+Wq2V9pOAz5ZVQcBn+ynN5JkO+As4GjgYOAVSQ4eWL4UeD5wyyz7IkmSNG/NNrQdB3ywf/9B4CUjyhwGrK6qm6vqYeCCvt6U9wBvBGqWfZEkSZq3Zhva9q2qdQD9v/uMKLMfsGZgem0/jyQvBr5ZVdfM1FCSk5KsTLJy/fr1s+y2JElSW7afqUCSfwKeOGLRWyZsIyPmVZKd+nW8YJKVVNU5wDkAy5cv96qcJElaUGYMbVX18+OWJflWkkVVtS7JIuD2EcXWAksHppcAtwJPBg4ErkkyNf9LSQ6rqts2YxskSZLmvdneHl0BvKp//yrgohFlvggclOTAJDsAxwMrquq6qtqnqpZV1TK6cHeIgU2SJGlTsw1t7wKen+RrdE+AvgsgyeIkFwNU1SPAKcClwCrgo1V1/SzblSRJWlBmvD06naq6EzhyxPxbgWMGpi8GLp5hXctm0xdJkqT5zL+IIEmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgMMbZIkSQ0wtEmSJDXA0CZJktQAQ5skSVIDDG2SJEkNMLRJkiQ1wNAmSZLUAEObJElSAwxtkiRJDTC0SZIkNSBVtbX7sNmS3A/cuLX7sQ3aC7hja3diG+OYjOa4jOa4jOa4bMoxGc1xGe1Hq2rX2a5k+7noyVZwY1Ut39qd2NYkWem4bMwxGc1xGc1xGc1x2ZRjMprjMlqSlXOxHm+PSpIkNcDQJkmS1IBWQ9s5W7sD2yjHZVOOyWiOy2iOy2iOy6Yck9Ecl9HmZFyafBBBkiRpoWn1SpskSdKCYmiTJElqwDYV2pIcleTGJKuTnDZieZKc2S+/Nskhk9Zt2QTjckI/Htcm+WySZwws+3qS65JcPVePHG8rJhiX5ya5t9/2q5P87qR1WzbBuPzWwJh8Jcn3k+zZL5uXx0uSc5PcnuQrY5Yv1HPLTOOy4M4tE4zJQj2vzDQuC+68ApBkaZJ/TrIqyfVJTh1RZu7OL1W1TbyA7YCbgCcBOwDXAAcPlTkGuAQIcDhw5aR1W31NOC4/A+zRvz96alz66a8De23t7dhK4/Jc4BNbUrfV1+ZuG3As8KkFcLw8BzgE+MqY5Qvu3DLhuCzEc8tMY7LgziuTjMtQ2QVxXum3bRFwSP9+V+CrP8zssi1daTsMWF1VN1fVw8AFwHFDZY4DPlSdzwO7J1k0Yd1WzbhtVfXZqrq7n/w8sORR7uPWMJt9vqCPlyGvAD7yqPRsK6qqzwB3TVNkIZ5bZhyXhXhumeBYGWdBHytDFsR5BaCq1lXVl/r39wOrgP2Gis3Z+WVbCm37AWsGptey6YaPKzNJ3VZt7rb9El2in1LAZUmuSnLSD6F/W8uk4/KsJNckuSTJ0zazbosm3rYkOwFHARcOzJ6vx8tMFuK5ZXMtlHPLJBbaeWViC/m8kmQZ8JPAlUOL5uz8si39GauMmDf8+0jGlZmkbqsm3rYkz6M7sf7swOxnV9WtSfYBLk9yQ/8TU+smGZcvAQdU1QNJjgH+Djhowrqt2pxtOxa4oqoGf3qer8fLTBbiuWViC+zcMpOFeF7ZHAvyvJJkF7qg+htVdd/w4hFVtuj8si1daVsLLB2YXgLcOmGZSeq2aqJtS/J04P3AcVV159T8qrq1//d24ON0l2PngxnHparuq6oH+vcXA49NstckdRu2Odt2PEO3MObx8TKThXhumcgCPLdMa4GeVzbHgjuvJHksXWA7v6o+NqLInJ1ftqXQ9kXgoCQHJtmBbsevGCqzAnhl/yTG4cC9VbVuwrqtmnHbkuwPfAw4saq+OjB/5yS7Tr0HXgCMfPKnQZOMyxOTpH9/GN3xfuckdRs20bYl2Q04ArhoYN58Pl5mshDPLTNaoOeWaS3Q88pEFuJ5pT8W/gpYVVXvHlNszs4v28zt0ap6JMkpwKV0T1ScW1XXJzm5X342cDHdUxirgYeA10xXdytsxpybcFx+F3gC8Of9ueSRqloO7At8vJ+3PfA3VfWPW2Ez5tyE4/ILwK8keQT4NnB8dY/sLPTjBeClwGVV9eBA9Xl7vCT5CN1Tf3slWQv8HvBYWLjnFphoXBbcuWWCMVlw5xWYaFxggZ1Xes8GTgSuS3J1P+/NwP4w9+cX/4yVJElSA7al26OSJEkaw9AmSZLUAEObJElSAwxtkiRJDTC0SZIkNcDQJkmS1ABDmyRJUgP+P4hebykYrAApAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "tags": [], "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "ouvumpEfjOCw" }, "source": [ "## Exact Solution\n", "The initial value problem has the exact solution\n", "$$y(t)=2e^{-t}+t-1.\\end{equation} \n", "The figure below plots the exact solution." ] }, { "cell_type": "code", "metadata": { "id": "g6Jc6lNxjOCw", "outputId": "adfdc4ab-cf46-482c-93f2-d64f5278e097" }, "source": [ "IC=1 # Intial condtion\n", "y=(IC+1)*np.exp(-t)+t-1\n", "fig = plt.figure(figsize=(6,4))\n", "plt.plot(t,y,'o-',color='black')\n", "plt.title('Exact Solution ')\n", "plt.xlabel('time')\n", "plt.show()" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmyUlEQVR4nO3de3SU9bX/8fdOCAgIRkoEAZPUSxXkV7Hipf7sUVfrqvW0SmtVEKgokoM3tAWpiFA4ipXaHrRKOUVUiqSiLeBS0V+laqUt3oKiUpVKkSBCCxYEFAyQ7N8fM8EhmVuSmXnm8nmtNcvMc90ZH3a+832+3/2YuyMiIrmvKOgAREQkNZTQRUTyhBK6iEieUEIXEckTSugiInlCCV1EJE8ooYu0gpmtM7NvtHLfcjP7xMyKUx2XFDYldMmIcALcHU5kja9703i+s8xsQ4Jt+pjZQjP7yMy2m9lbZjYiDbEckPzdfb27H+zu9ak+lxS2dkEHIAXlO+7+x6CDiPAQ8AZQAdQB/wfoGWhEIm2gFroEzsxmmdnvI95PN7NnLeRQM3vSzLaY2bbwz30itu1mZg+a2cbw+sfMrDPwNNAr4ttAryinPhmY6+6fuvs+d3/d3Z+OOPb5ZvY3M/vYzP5kZn1jxD/XzG6LeL//24GZPQSUA0+E4xhvZpVm5mbWLrxNLzN73My2mtkaMxsVcawpZvaomc0zs53heAa29rOW/KaELtlgLPBlMxthZl8DRgKXeaguRRHwIKFWdDmwG4jsqnkI6AQcDxwGzHD3T4FvARvDXRsHu/vGKOd9CZhpZoPNrDxyhZl9CXgYuAEoA54ilJTbt+QXc/fhwHpC304OdvefRdnsYWAD0Av4PnC7mX09Yv35wAKgFHi8ye8vsp8SumTSY+HWbuNrFIC77wKGAf8DzAeuc/cN4XX/dveF7r7L3XcC04AzAczscEKJe7S7b3P3ve7+QgviuQj4MzAJeN/MVprZyeF1lwBL3H2pu+8Ffg50BE5v42dwADM7AjgD+LG7f+buK4E5wPCIzf7i7k+F+9wfAk5IZQySP5TQJZMGuXtpxOu+xhXu/gqwFjDg0cblZtbJzH5tZrVmtgNYBpSGR4gcAWx1922tCSb8R+Amdz8e6AGsJPRHxwi1lmsjtm0APgB6t+ZccfQi9DvsjFhW2+Q8/4z4eRdwUGN3jUgkJXTJCmZ2DdAB2AiMj1g1FjgWONXduwL/0bgLoQTbzcxKoxyyRWVE3f0jQq3wXkC3cBwVEfEZoT8gH0bZ/VNC3T6Nmt5YjRfLRkK/Q5eIZeUxziMSlxK6BC7cX30boW6X4cB4MxsQXt2FUL/5x2bWDfhJ437uvonQzc9fhW+elphZY8L/F/AFMzskznmnm1l/M2sXTqhXAWvc/d+EviX8p5l93cxKCP1hqQOWRznUSuC88A3anoT63SP9CzgyWgzu/kH4mD81s4PM7MuE7iFUx4pbJBYldMmkxpEeja/F4a6D+cB0d3/D3d8DbgYeMrMOwF2E+q4/InQT8/81OeZwYC/wLrCZcDJ193cJ3WxcG+6vjzbKpROwGPiYUHdPBaEbkLj7akJ/YO4Jn/s7hG5s7olynMbhj+uAZ4BHmqz/KXBLOI5xUfYfAlQSaq0vBn7i7kujbCcSl+kBFyIi+UEtdBGRPKGELiKSJ5TQRUTyhBK6iEieCGxyQvfu3b2ysjKo04uI5KQVK1Z85O5l0dYFltArKyupqakJ6vQiIjnJzGpjrVOXi4hInlBCFxHJE0roIiJ5QgldRCRPKKGLiOQJJXQRkQyprq6msrKSoqIiKisrqa5ObVFNFckXEcmA6upqqqqq2LVrFwC1tbVUVVUBMHTo0JScQy10EZEMmDhx4v5k3mjXrl1MnDgxZedQQhcRyYD169e3aHlrKKGLiGRAnz59oi4vLy9P2TmU0EVEMuDkk09utqxTp05MmzYtZedQQhcRSbM1a9awZMkSTjvtNCoqKjAzKioqmD17dspuiIJGuYiIpJW7c80119ChQwcWLlxIr17RHm+bGkroIiJp9Lvf/Y5nnnmGX/7yl2lN5hDgQ6IHDhzoKp8rIvls+/bt9O3bl169evHyyy9TXFzc5mOa2Qp3HxhtnVroIiJpMmnSJP75z3/y+OOPpySZJ6KboiIiabBixQpmzpzJ1VdfzcCBURvUKaeELiKSYvX19YwePZqysjJuu+22jJ1XCV1EJEUai2+1a9eOmpoaLrzwQkpLSzN2fiV0EZEUaCy+VVv7+SM/586dm/KKivEooYuIpEAmim8looQuIpICmSi+lYgSuohICmSi+FYiSugiIilwwgknNFuW6uJbiSihi4i00ZtvvsnTTz/NmWeemdbiW4lopqiISBvU19dTVVVFt27dWLRoEd26dQssFiV0EZE2+PWvf83LL7/M/PnzA03moC4XEZFW27hxIxMmTOCcc87h0ksvDTocJXQRkdYaM2YMe/bsYdasWZhZ0OEkTuhm9oCZbTazVTHWDzWzN8Ov5WbW/FaviEieeeKJJ1i4cCGTJ0/mqKOOCjocILkW+lzg3Djr3wfOdPcvA7cCs1MQl4hIVqqurqa8vJzzzz+fkpISevfuHXRI+yW8Keruy8ysMs765RFvXwKij64XEclxjfVaGqf47927l6uuuori4uKMDk+MJdV96COBp2OtNLMqM6sxs5otW7ak+NQiIumVDfVa4klZQjezswkl9B/H2sbdZ7v7QHcfWFZWlqpTi4hkRDbUa4knJQndzL4MzAEucPd/p+KYIiLZ5pBDDom6PJP1WuJpc0I3s3JgETDc3f/e9pBERLLPe++9xyeffNLs2aCZrtcSTzLDFh8GXgSONbMNZjbSzEab2ejwJpOBLwC/MrOVZlaTxnhFRDKuoaGBUaNG0blzZ+6+++5A67XEY+4eyIkHDhzoNTXK/SKS/ebMmcOoUaO47777uPLKKwONxcxWuHvUp05rpqiISBybNm1i3LhxnHXWWYwcOTLocOJSQhcRieO6666jrq6O2bNnZ8X0/nhUbVFEJIbFixezcOFC7rjjDo455pigw0lIfegiIlF8/PHH9OvXjx49evDKK69QUlISdEiA+tBFRJJWXV1NZWUlhx56KJs2beLCCy/MmmSeiBK6iEhYY62W2tra/ct++tOfUl1dHWBUyVNCFxEJy/ZaLYkooYuIhGV7rZZElNBFRMJ69OgRdXm21GpJRAldRASoq6ujuLi42VjzbKrVkogSuogIMG3aND788EPGjh2btbVaEtE4dBEpeCtXruTkk0/m0ksv5Te/+U3Q4cSlcegiIjHs3buXK664gi984QvMmDEj6HDaRFP/RaSg3Xnnnbz++ussWrSIbt26BR1Om6iFLiIF6+2332bq1KlcfPHFfPe73w06nDZTQheRgtI4tb+oqIgTTzyRkpIS7rnnnqDDSgkldBEpGJFT+92dPXv2sHfvXpYuXRp0aCmhhC4iBSPa1P49e/bkzNT+RJTQRaRg5PrU/kSU0EWkYMSawp8rU/sTUUIXkYJx7bXXNluWS1P7E1FCF5GCUF9fz8KFC+ncuTN9+vTJyan9iWhikYgUhF/84he89NJL/Pa3v2XIkCFBh5MWaqGLSN57++23mTRpEt/73vcYPHhw0OGkjRK6iOS1ffv2cdlll9G1a1dmzZrVrDxuPlGXi4jktZ/97GfU1NTw6KOPcthhhwUdTlolbKGb2QNmttnMVsVYf5yZvWhmdWY2LvUhioi0zptvvsmUKVO45JJLuOiii4IOJ+2S6XKZC5wbZ/1WYAzw81QEJCLSVtXV1VRUVHDCCSdQX1/PWWedFXRIGZEwobv7MkJJO9b6ze7+KrA3lYGJiLRGY72WxtmfDQ0NjB07lurq6oAjS7+M3hQ1syozqzGzmi1btmTy1CJSIKLVa9m1a1fe1GuJJ6MJ3d1nu/tAdx9YVlaWyVOLSIHI93ot8WjYoojklYMPPjjq8nyp1xKPErqI5I3nnnuOnTt30q7dgSOy86leSzzJDFt8GHgRONbMNpjZSDMbbWajw+t7mtkG4EfALeFtuqY3bBGRA3388ceMGDGCL33pS8yePZuKioq8rNcST8KJRe4et+iBu/8T6JOyiEREWmHMmDFs3LiR5cuXc8opp3D55ZcHHVLGqctFRHLewoULeeihh5g4cSKnnHJK0OEERgldRHLapk2b+K//+i9OOukkbrnllqDDCZQSuojkLHfnyiuv5NNPP+Whhx6ipKQk6JACpYQuIjmnurqayspKioqKeOqpp/j+979P3759gw4rcEroIpJTGqf219bW7l+2aNGigpjan4gSuojklEKe2p+IErqI5JRCntqfiBK6iOSUWA+pKISp/YkooYtIztixYwf19fXNHiNXKFP7E1FCF5Gccd1117F161YmTZpUkFP7E9EzRUUkJyxYsIB58+YxefJkpk6dytSpU4MOKeuohS4iWW/9+vWMHj2a0047jUmTJgUdTtZSQheRrFZfX8/w4cOpr69n/vz5zUrjyueU0EUkKzXOBm3Xrh3Lli3j0ksv5aijjgo6rKymhC4iWSfabND58+drNmgCSugiknU0G7R1lNBFJOtoNmjrKKGLSNbp1q1b1OWaDRqfErqIZJU1a9bwySefUFR0YHrSbNDElNBFJGvU1dVxySWX0KlTJ2bMmKHZoC2kAZ0ikjUmTJjAa6+9xuLFixk0aBBjxowJOqScoha6iGSFJUuWMGPGDK699loGDRoUdDg5SQldRAL34YcfMmLECE444QTuvPPOoMPJWUroIhKo+vp6hg0bxq5du3jkkUc46KCDgg4pZymhi0ggIqf2/+lPf2LYsGEce+yxQYeV0xImdDN7wMw2m9mqGOvNzH5pZmvM7E0z+0rqwwyJfNJ3ZWWlpgGL5ChN7U+PZFroc4Fz46z/FnBM+FUFzGp7WM1FXgDuTm1tLVVVVboARHKQpvanR8KE7u7LgK1xNrkAmOchLwGlZnZ4qgJspAtAJH9oan96pKIPvTfwQcT7DeFlzZhZlZnVmFnNli1bWnQSXQAi+aNr165Rl2tqf9ukIqFblGUebUN3n+3uA919YFlZWYtOEut/tC4Akdzy/PPPs337doqLiw9Yrqn9bZeKhL4BOCLifR9gYwqOe4Bp06bRqVOnA5aVlJToAhDJIZs2bWLIkCEcd9xxzJ49W1P7UywVU/8fB641swXAqcB2d9+UguMeoPF/9MSJE1m/fj0dOnSgffv2XHjhhak+lYikwb59+xgyZAg7duzgj3/8I/379+eKK64IOqy8ksywxYeBF4FjzWyDmY00s9FmNjq8yVPAWmANcB9wdbqCHTp0KOvWraOhoYElS5awY8cO5s2bl67TiUgKTZkyhRdeeIFZs2bRv3//oMPJS+Yetbs77QYOHOg1NTWt3t/dOfXUU9m6dSurV69u1h8nItnj6aef5rzzzmPkyJHMmTMn6HBympmtcPeB0dbl7ExRM+Omm27iH//4BwsXLgw6HBFpInIi4Le//W2OOOII7rnnnqDDyms5m9ABLrjgAr70pS9xxx13ENQ3DRFprulEwIaGBrZs2cKiRYuCDi2v5XRCLy4uZvz48bz++uv88Y9/DDocEQmLNhHws88+00TANMvZPvRGdXV1HHnkkRx33HE8++yzKYhMRNqqqKgo6rdmM6OhoSGAiPJHXvahN+rQoQM/+tGPeO6553j11VeDDkdEgJ49e0ZdromA6ZXzCR2gqqqK0tJSpk+fHnQoIgVv27Zt1NfXY3bgJHLNBE2/vEjoXbp04ZprrmHRokWsXr066HBEClZDQwPDhg1j27ZtTJ48WTNBMywvEjrAmDFjKC4u5qSTTlK9dJGA3HrrrTz11FPcddddTJkyZf9EwHXr1imZZ0Aqpv5nhaVLlwLw6aefAuyvlw7oQhLJgCeffJIpU6bwgx/8gKuuuirocApSzo9yaVRZWXnA008aVVRUsG7dupSdR0Q+V11dvb++EsARRxzBu+++S8eOHQOOLH/l9SiXRqqXLpJZTScPubsmDwUsbxK66qWLZFa0yUO7d+/W5KEA5U1Cj1YvvWPHjhomJZIm+lacffImoQ8dOvSAgvkAp556qm6IiqRJ9+7doy7Xt+Lg5E1ChwPrpY8ZM4Y///nPvPfee0GHJZJ33nnnHXbu3KnJQ1kmrxJ6pJtvvpkOHTowZcqUoEMRySvbtm3j/PPPp0uXLtx1112aPJRF8jah9+jRgzFjxvDwww+zatWqoMMRyQv79u3jkksuoba2lkWLFjFmzBhNHsoieZvQAW688Ua6dOnCpEmTgg5FJC+MHz+epUuXMmvWLM4444ygw5Em8jqhd+vWjXHjxvHYY4+pEqNIGz344IPMmDGDMWPGMHLkyKDDkSjyOqED3HDDDXTv3p1bbrkl6FBEck7kY+SuuOIKjj/+eH7xi18EHZbEkPcJvUuXLtx0000888wzLFu2LOhwRHJG05mgAGvXruWRRx4JODKJJW9qucSze/dujj76aI488kiWLVvWbKiViDSn+kjZqSBqucTTsWNHbrnlFv7yl7/Qs2dPldcVSYJmguaegkjoEJrwYGZs3rwZd99fXldJXSS6gw8+OOpyzQTNXgWT0H/yk580e2jtrl27VEhIJIp7772XnTt30q7dgY9M0EzQ7JZUQjezc81stZmtMbOboqw/1MwWm9mbZvaKmfVPfahto6+PIslZsmQJ119/Peeffz4PPPCAZoLmkIRPLDKzYmAmcA6wAXjVzB5397cjNrsZWOnu3zWz48Lbfz0dAbdWeXl51Bs8+voo8rk33niDwYMHM2DAAH7729/SuXNnhg8fHnRYkqRkWuinAGvcfa277wEWABc02aYf8CyAu78LVJpZj5RG2kYqryvSXOQ48z59+nD22WdTWlrKE088QefOnYMOT1oomYTeG/gg4v2G8LJIbwDfAzCzU4AKoE/TA5lZlZnVmFnNli1bWhdxK0UrrztgwAB9fZSC1XSc+Ycffsi2bdu4+uqr6dWrV9DhSSskk9CjDdpuOnj9DuBQM1sJXAe8DuxrtpP7bHcf6O4Dy8rKWhprm0WW1500aRIvvvgif/3rXzMeh0g2iPbEIYBf//rXAUQjqZBMQt8AHBHxvg+wMXIDd9/h7pe7+wDgB0AZ8H6qgkyHH//4x/Tp04cxY8ZQX18fdDgiGaeBAvknmYT+KnCMmX3RzNoDg4HHIzcws9LwOoArgWXuviO1oaZW586d+dnPfsZrr73Ggw8+GHQ4Ihmn5/Dmn4QJ3d33AdcCfwDeAR5197+Z2WgzGx3erC/wNzN7F/gWcH26Ak6lwYMHc8YZZ3DzzTezffv2oMMRyahzzjmn2TKNM89tBVHLJZ7XXnuNgQMH8sMf/lBV5KRgPPnkkwwaNIh+/fqxfft2PvjgA8rLy5k2bZoGCmS5eLVcCj6hA4waNYq5c+fy1ltvcdxxxwUdjkhavfjii3z961/n+OOP5/nnn485xV+yU8EX50pk2rRplJSUcOKJJ6pwl+S1d999l29/+9v07t2bJUuWKJnnGSV0YOnSpezbt4/PPvtMhbsk70ROHurfvz/79u3jD3/4A4cddljQoUmKKaETGo+7d+/eA5apcJfkg6aTh+rr69mzZw8vvvhi0KFJGqgPHSgqKmpWiRHAzGhoaAggIpHU0EMq8o/60BPQeFzJV5o8VFiU0IleuMvMmDp1akARibTdvn37OOigg6KuU2MlPymh07xwV1lZ2f6boyK5qKGhgZEjR7J7925KSkoOWKfJQ/lLCT0ssnDX5s2bufTSS7nttttYtWpV0KGJtIi7c9111zFv3jz++7//mwcffFAPqSgQuikaw0cffUS/fv2orKxk+fLlzR7FJZKN3J0JEyYwffp0brzxRqZPn76/XLTkB90UbYXu3btz77338uqrrzJjxoygwxGJKXKc+aGHHsr06dO56qqrlMwLkBJ6HBdddBGDBg1i8uTJ/P3vfw86HJFmmo4z3759O8XFxZx++ulK5gVIXS4JbNq0iX79+tG/f39eeOEFior0N1Cyh8aZFx51ubTB4Ycfzl133cVf/vIXunfvrlovklU0zlwiKaEnoV27dhQVFbFt2zbVepGscuihh0ZdrnHmhUkJPQkTJ05sVgJAtV4kaDNnzmTr1q0UFxcfsFzjzAuXEnoS9LVWss0999zDtddeywUXXMD999+vceYCgAZXJ6G8vDzqjSd9rZVMqa6uZuLEiaxfv57S0lK2bdvGoEGDeOSRR2jfvj2XXXZZ0CFKFlALPQnRar0AjB8/PoBopNA0HZq4bds2iouL+e53v0v79u0TH0AKhhJ6EprWejn88MMpKSlh4cKF1NfXBx2e5LmJEyeya9euA5bV19czefLkgCKSbKWEnqTIWi8bN27kf//3f3nuuee4/fbbgw5N8pzu4UiylNBb6fLLL2fYsGFMmTKFF154IehwJE81NDTEfO6n7uFIU0rorWRmzJo1i6OPPpohQ4awefPmoEOSPLN3715GjBjBzp07mxWH09BEiUYJvQ0OPvhgHn30UTZv3kxFRYVmkUrK7N69mwsvvJCHHnqIW2+9lblz52pooiSkYYtttGrVKoqLi/nss88A9s8iBfQPTlplx44dnH/++SxbtoyZM2dy9dVXA7qeJLGkWuhmdq6ZrTazNWZ2U5T1h5jZE2b2hpn9zcwuT32o2WnixIns2bPngGWaRSotFVkCt3v37ixbtozq6ur9yVwkGQlb6GZWDMwEzgE2AK+a2ePu/nbEZtcAb7v7d8ysDFhtZtXuvifKIfOKRiBIWzWOM28cmrh37146dOjQrNyESCLJtNBPAda4+9pwgl4AXNBkGwe6WKgA88HAVmBfSiPNUrFGGvTu3TvDkUiuijbOvK6uTt/ypMWSSei9gQ8i3m8IL4t0L9AX2Ai8BVzv7s2aF2ZWZWY1ZlazZcuWVoacXWLNIu3QoQO7d+8OICLJNfqWJ6mSTEKP9tiTpk/F+CawEugFDADuNbOuzXZyn+3uA919YFlZWQtDzU5NZ5FWVFRw/fXXs3btWkaMGKGvzRLXvffeS6yHzGicubRUMqNcNgBHRLzvQ6glHuly4A4PXZlrzOx94DjglZREmeWGDh3abARC7969GT9+PEcffbTGC0sz9fX1jB07lrvvvpuvfOUrvPPOOwd8o9M4c2mNZFrorwLHmNkXzaw9MBh4vMk264GvA5hZD+BYYG0qA80148aNY9SoUdx+++08+OCDQYcjAYscxVJeXs4pp5zC3XffzQ9/+ENeeeUV7rvvPo0zl7Zz94Qv4Dzg78A/gInhZaOB0eGfewHPEOo/XwUMS3TMk046yfPdnj17/JxzznEz88MOO8zNzCsqKnz+/PlBhyYZNH/+fO/UqZMT6qrc/7rsssuCDk1yEFDjMfKqHhKdZnPmzKGqquqAftJOnTqpBVZA9CBnSaV4D4lWQk8z/WOWoqKiqDc+zUw3zaXF4iV01XJJMw1JK2x79uyhc+fOUddpFIukmhJ6msX6R9urV68MRyKZtnHjRs4++2w++eQTVUuUjFBCT7NYE4/q6urU5ZJnIkey9OzZk759+7Jy5UoWLFigaomSGbHulqb7VQijXBrNnz/fKyoq9o9yufXWW720tNQrKyv9/fffDzo8SYFoI1nMzO+4446gQ5M8g0a5ZJ8VK1bwjW98g9LSUm644QZmzJjB+vXrKS8vZ9q0aWq95Rjd/JZM0SiXLLVixQq+9rWv8dlnn2lYY47TSBbJFI1yyVInnXQShxxySLNEoHrquaO+vp5p06apHotkBSX0gP3rX/+KulzDGrNfbW0tZ599NrfccgunnXYaHTt2PGC9RrJIpimhByxWC04tu+wTOYqlrKxs/yiWefPmsXz5ctVjkcApoQcs1rDGE088UX2vWaTxqUK1tbW4Ox999BF1dXVMnTqV4cOHY2YMHTqUdevW0dDQwLp165TMJeOU0APWtJ56eXk5Z555Jo899hiDBw9u9iQbCUa0pwo1NDRw9913BxSRSHNK6FkgsmVXW1vL888/z5133snvf/97+vfvT58+fSgqKqKyspLq6uqgwy04mzZtijokEXSvQ7KLEnoWMjPGjRvHDTfcwPvvv8+HH36Iu1NbW0tVVZWSehpF9pNXVFRw5ZVX0rdv35jb616HZBMl9Cy2aNGiZss0pDF9mvaTr1+/nvvvv59evXrx85//vNm9Do1ikWyjhJ7FVKkxs6L1kwN8+umnjB07ttmzYzWKRbKNZopmsVjTyYuKiliwYAEXXXRRAFHlp4aGBoqLi6Ou02xPySaaKZqjog1pPOigg6ioqODiiy9m2LBh3Hffffv7fHXTtHWee+45Tj755Jjr1U8uOSNW1a50vwqp2mJbNK3UOH/+fN+zZ49PmTLFzczN7IAKf506ddIzS+OI/DwPP/xwHzBggANeXl7uV111VbOKifo8JdsQp9qiEnoO69mzZ7MHDwNeUVERdGhZKdbDmgcPHuy7d+/ev03TP6Ai2SReQlcfeg5Thb+W6d27Nxs3bmy2XCVuJZeoDz1PxerbNTNmz57Nvn37DhhXXQh97NF+3xUrVnDxxRdHTeagUUOSR2I13dP9UpdL20XrQjjooIP8mGOOccB79erlHTp0KJg+4WifR1FRkQPetWtX79q1q7qoJOcRp8tFLfQc1rQOTEVFBXPmzGH16tUsXryYzZs3U1dXd8A++TwxKVa9ldLSUj744AN+9atfaXKQ5LdYmT7dL7XQ06/pCJjGl5m5e+7eAGwa97x58/yZZ56J+rtG/r7R9s2V31mkEW0d5QKcC6wG1gA3RVl/I7Ay/FoF1APd4h1TCT39KioqYia5r371qznZHRPrYcyR/236UpeK5JN4CT1hl4uZFQMzgW8B/YAhZtavSSv/Tncf4O4DgAnAC+6+NbnvCJIusSYmffOb3+Sll15K2B0T1A3VWOfdt28fY8eObdat4u50796dBx54QF0qUthiZXr/vPX9VeAPEe8nABPibP9bYFSi46qFnhmxuhhitWYBf+KJJ/z+++9P6ySbWHFFa4G3b9/eTz75ZC8tLU3YraIuFcl3tKXLBfg+MCfi/XDg3hjbdgK2kqC7xZXQAxerOyZeoiei+yJR4oy3PlrS7tSpk8+aNcsPO+ywqOctLi72K664wrt3765uFSlobU3oF0VJ6PfE2PYS4Ik4x6oCaoCa8vLyDP36Ek2spDp37lx/9tln4yb1M88800tKSg5Y1rFjR//Nb34T99gzZ870d955x3v06BH3+Ila4JqeL4WsrQk96S4XYDFwaaJjulroWSFeKzpWC75Dhw77x3ZHe7Vr1y5hKz/eK1YLPbIFrm4VKWRtTejtgLXAF4H2wBvA8VG2OyTc3dI50TFdCT3rxWsJx0vYEyZMiJuwH3744bhJWy1wkfjalNBD+3Me8HfgH8DE8LLRwOiIbUYAC5I5niuh54RYLeFYrffGVnSi9YmStlrgIrG1OaGn46WEnruSSciJWtlK2iKto4QuKdeWUS4i0nrxErrK54qI5BCVzxURKQBK6CIieUIJXUQkTyihi4jkCSV0EZE8EdgoFzPbAtS2cvfuwEcpDCdVsjUuyN7YFFfLKK6Wyce4Kty9LNqKwBJ6W5hZTaxhO0HK1rgge2NTXC2juFqm0OJSl4uISJ5QQhcRyRO5mtBnBx1ADNkaF2RvbIqrZRRXyxRUXDnZhy4iIs3lagtdRESaUEIXEckTWZfQzexcM1ttZmvM7KYo683Mfhle/6aZfSXZfdMc19BwPG+a2XIzOyFi3Toze8vMVppZSktMJhHXWWa2PXzulWY2Odl90xzXjRExrTKzejPrFl6Xzs/rATPbbGarYqwP6vpKFFdQ11eiuIK6vhLFlfHry8yOMLPnzewdM/ubmV0fZZv0Xl+x6uoG8QKKCT0V6Ug+f9xdvybbnAc8DRhwGvBysvumOa7TgUPDP3+rMa7w+3VA94A+r7OAJ1uzbzrjarL9d4Dn0v15hY/9H8BXgFUx1mf8+koyroxfX0nGlfHrK5m4gri+gMOBr4R/7kLoKW8ZzV/Z1kI/BVjj7mvdfQ+wALigyTYXAPM85CWg1MwOT3LftMXl7svdfVv47UtAnxSdu01xpWnfVB97CPBwis4dl7svI/Ts21iCuL4SxhXQ9ZXM5xVLoJ9XExm5vtx9k7u/Fv55J/AO0LvJZmm9vrItofcGPoh4v4HmH0isbZLZN51xRRpJ6K9wIweeMbMVZlaVophaEtdXzewNM3vazI5v4b7pjAsz6wScCyyMWJyuzysZQVxfLZWp6ytZmb6+khbU9WVmlcCJwMtNVqX1+mrX0h3SzKIsazquMtY2yezbWkkf28zOJvQP7oyIxf/X3Tea2WHAUjN7N9zCyERcrxGq/fCJmZ0HPAYck+S+6Yyr0XeAv7p7ZGsrXZ9XMoK4vpKW4esrGUFcXy2R8evLzA4m9AfkBnff0XR1lF1Sdn1lWwt9A3BExPs+wMYkt0lm33TGhZl9GZgDXODu/25c7u4bw//dDCwm9PUqI3G5+w53/yT881NAiZl1T2bfdMYVYTBNvg6n8fNKRhDXV1ICuL4SCuj6aomMXl9mVkIomVe7+6Iom6T3+kr1jYG2vAh9Y1gLfJHPbwwc32Sb/+TAmwqvJLtvmuMqB9YApzdZ3hnoEvHzcuDcDMbVk88nkJ0CrA9/doF+XuHtDiHUD9o5E59XxDkqiX2TL+PXV5JxZfz6SjKujF9fycQVxPUV/r3nAXfF2Sat11fKPtwU/k86j9Dd4X8AE8PLRgOjIz60meH1bwED4+2bwbjmANuAleFXTXj5keH/OW8AfwsgrmvD532D0M200+Ptm6m4wu9HAAua7Jfuz+thYBOwl1CraGSWXF+J4grq+koUV1DXV9y4gri+CHWDOfBmxP+n8zJ5fWnqv4hInsi2PnQREWklJXQRkTyhhC4ikieU0EVE8oQSuohInlBCl4JhZqVmdnX4515m9vugYxJJJQ1blIIRrq/xpLv3DzoWkXTItlouIul0B3CUma0E3gP6unt/MxsBDCJUwrQ/8AtCs/WGA3XAee6+1cyOIjQppAzYBYxy93cz/UuIxKIuFykkNwH/cPcBwI1N1vUHLiU0fX0asMvdTwReBH4Q3mY2cJ27nwSMA36ViaBFkqUWukjI8x6qYb3TzLYDT4SXvwV8OVxB73Tgd2b7C+N1yHyYIrEpoYuE1EX83BDxvoHQv5Mi4ONw614kK6nLRQrJTkKPBmsxD9W1ft/MLoL9z4Y8IZXBibSVEroUDA/VEP9r+MHCd7biEEOBkWbWWKkvZY9UE0kFDVsUEckTaqGLiOQJJXQRkTyhhC4ikieU0EVE8oQSuohInlBCFxHJE0roIiJ54v8D2VbzXO4HRHIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "tags": [], "needs_background": "light" } } ] }, { "cell_type": "code", "metadata": { "id": "bWlRDTDnjOCz" }, "source": [ "# Initial Condition\n", "w=np.zeros(N+1)\n", "#np.zeros(N+1)\n", "w[0]=IC" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "O6yWJtXLjOC1" }, "source": [ "\n", "\n", "## 2-step Abysmal Butler Method \n", "\n", "The 2-step Abysmal Butler difference equation is\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{2}(4f(t_i,w_i)-3f(t_{i-1},w_{i-1})).\\end{equation} \n", "\n", "For $i=0$ the system of difference equation is:\n", "\\begin{equation}w_{1} = w_{0} + \\frac{h}{2}(4(t_0-w_0)-3(t_{-1}-w_{-1})) \\end{equation} \n", "this is not solvable as $w_{-1}$ is unknown.\n", "\n", "For $i=1$ the difference equation is:\n", "\\begin{equation}w_{2} = w_{1} + \\frac{h}{2}(4(t_1-w_1)-3(t_{0}-w_{0})),\\end{equation}\n", "this is not solvable as $w_{1}$ is unknown. $w_1$ can be approximated using a one step method. Here, as the exact solution is known,\n", "\\begin{equation}w_1=2e^{-t_1}+t_1-1.\\end{equation} \n" ] }, { "cell_type": "code", "metadata": { "id": "AbLbG_1AjOC2" }, "source": [ "### Initial conditions\n", "w=np.zeros(len(t))\n", "w0=np.zeros(len(t))\n", "w[0]=IC\n", "w[1]=y[1]" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "AMuQXLLXjOC5" }, "source": [ "### Loop" ] }, { "cell_type": "code", "metadata": { "id": "_lBut4d0jOC5" }, "source": [ "for k in range (1,N):\n", " w[k+1]=(w[k]+h/2.0*(4*myfun_ty(t[k],w[k])-3*myfun_ty(t[k-1],w[k-1]))) \n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "ANL1cpKmjOC7" }, "source": [ "### Plotting solution" ] }, { "cell_type": "code", "metadata": { "id": "n1_2486MjOC8" }, "source": [ "def plotting(t,w,y):\n", " fig = plt.figure(figsize=(10,4))\n", " plt.plot(t,y, 'o-',color='black',label='Exact')\n", " plt.plot(t,w,'^:',color='red',label='Abysmal-Butler')\n", " plt.xlabel('time')\n", " plt.legend()\n", " plt.show " ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Vr__klcUjOC-" }, "source": [ "The plot below shows the exact solution (black) and the Abysmal-Butler approximation (red) of the intial value problem.\n", "\n", "The Numerical approximation does not do a good job of approximating the exact solution and that is because it is inconsistent." ] }, { "cell_type": "code", "metadata": { "id": "rKNNxgLDjOC-", "outputId": "e39c0477-b50d-4eba-db8c-3d898bf8996d" }, "source": [ "plotting(t,w,y)" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEGCAYAAABB6hAxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABG0UlEQVR4nO3de5xN9f7H8dd3zF1CCDFmdHE5cVyaFCmphHA6dApNF1TSkeiUUhMK0905p06ikduPIRFdpVK6kArlFqVchmHc74a57e/vj2Wbiz33y94z834+HvOYvdZ3rb2+a/asvT/7u77fz9dYaxERERGRwvHzdgVEREREyjIFUyIiIiJFoGBKREREpAgUTImIiIgUgYIpERERkSLw99aBa9asaSMiIrx1eBEREZF8W7169QFrbS1PZV4LpiIiIli1apW3Di8iIiKSb8aY+JzKdJtPREREpAgUTImIiIgUgYIpERERkSLwWp8pT1JTU0lISOD06dPerooUQnBwMPXr1ycgIMDbVRERESk1PhVMJSQkUKVKFSIiIjDGeLs6UgDWWg4ePEhCQgINGzb0dnVERERKjU/d5jt9+jQ1atRQIFUGGWOoUaOGWhVFRKTUxMXFERERgZ+fHxEREcTFxXmlHj7VMgUokCrD9NqJiEhpiYuLY+DAgSQlJQEQHx/PwIEDAYiKiirVuvhUy5SIiIhIfkRHR58NpNySkpKIjo4u9boomMqmUqVKtGzZ8uzPiy++WGzPvWbNGhYtWlRszyciIlJR7dixo0DrS1KZDqZK4l5pSEgIa9asOfszYsSIYqipQ8GUiIhI0VlrqVKliseyBg0alHJtynAw5b5XGh8fj7X27L3Skuh8dvToURo3bszvv/8OQN++fZk8eTIADz30EJGRkVx++eWMHj367D4rV66kXbt2tGjRgjZt2nD06FFGjRrF3LlzadmyJXPnzi32eoqIiJR31loee+wxjh07hr9/1q7foaGhxMTElHqdfK4DutuwYcNYs2ZNjuU//PADycnJWdYlJSVx3333nQ10smvZsiX//e9/cz3uqVOnaNmy5dnlp556it69e/PGG2/Qr18/hg4dyuHDh3nggQcAiImJ4YILLiA9PZ0bb7yRdevW0aRJE3r37s3cuXO58sorOXbsGKGhoYwZM4ZVq1bxxhtv5OtvICIiIhlcLhePPPIIEyZM4JFHHqFNmzZER0ezY8cOGjRoQExMTKl3PgcfDqbykj2Qymt9frlv82XXqVMn5s2bx+DBg1m7du3Z9e+++y6xsbGkpaWRmJjIxo0bMcZQt25drrzySgDOP//8ItVJRESkonO5XAwaNIjJkyfz+OOP8/LLL2OM8UrwlJ3PBlN5tSBFREQQH3/uBM7h4eF8/fXXxV4fl8vFpk2bCAkJ4dChQ9SvX59t27bx6quvsnLlSqpXr06/fv04ffo01lqlCRARESkm6enp3H///UyfPp2nn36acePG+dTnbJntMxUTE0NoaGiWdSV5r/Q///kPTZs2Zc6cOQwYMIDU1FSOHTtG5cqVqVq1Knv37uXTTz8FoEmTJuzevZuVK1cCcPz4cdLS0qhSpQrHjx8vkfqJiIiUR2lpadxzzz1Mnz6d5557zucCKSjDwVRUVBSxsbGEh4djjCE8PJzY2NgiN/e5+0y5f0aMGMHmzZt5++23GT9+PNdeey3XXXcd48aNo0WLFrRq1YrLL7+cAQMGcM011wAQGBjI3LlzGTJkCC1atKBTp06cPn2ajh07snHjRnVAFxERyYfU1FTuvPNOZs+ezQsvvMCoUaN8LpACMNZarxw4MjLSrlq1Ksu6TZs20bRpU6/UR4qHXkMRESkOycnJ9O7dmw8++IDx48fzr3/9y6v1McasttZGeirz2T5TIiIiUjGdPn2a2267jUWLFvG///2Phx9+2NtVylWZvc0nIiIi5UfmRNzVq1dn0aJFvPXWWz4fSIFapkRERMTLsk9afPr0aQIDA6lcubKXa5Y/apkSERERr/I0aXFKSopXJi0uDAVTIiIi4lW+NGlxYSiYEhEREa/Zs2fPOXPsuXlj0uLCUDDlwcKFCzHG8NtvvwHw9ddf0717dy/XCs477zyP6yMiImjevDktW7akefPmfPDBB3k+1/Tp09m9e/fZ5euvv57sqSpERERK0s6dO+nQoQMAQUFBWcq8NWlxYZT9YCoxETp0gD17iu0p58yZQ/v27XnnnXeK7TlL2tKlS1mzZg3z58/nkUceyXP77MFUQaWnpxd6XxERkS1btnDttdeyZ88eli5dypQpU4o9EXdpKfvB1NixsGyZ87sYnDhxguXLlzNlypQswdSxY8fo2bMnf/nLXxg0aBAul4spU6bw6KOPnt1m8uTJ/Otf/+LkyZN069aNFi1a0KxZs7PZziMiInj66adp27YtkZGR/Pzzz3Tu3JlLLrmESZMmnT3+jTfeSOvWrfPdypTZsWPHqF69OgDbt2+nWbNmZ8teffVVnn32WebPn8+qVauIioqiZcuWnDp1KstzfP7557Rt25bWrVtz++23c+LEibP1HzNmDO3bt2fevHkFqpeIiIjbpk2buO666zh+/DhfffUV11xzDVFRUWzfvh2Xy8X27dvLTCAFvh5MXX89TJ/uPE5NdZZnzXKWk5KgbVt4+21wuWDqVGjXDhYscMoPHHC2/+gjZzmfLVfvv/8+Xbp0oVGjRlxwwQX8/PPPAPz000+MHz+e9evXs2XLFhYsWECfPn348MMPSU1NBWDatGn079+fxYsXc9FFF7F27Vo2bNhAly5dzj5/WFgYK1as4Nprr6Vfv37Mnz+fH374gVGjRgEQHBzMwoUL+fnnn1m6dCmPPfYY+clS37FjR5o1a0aHDh0YN25crtv+4x//IDIykri4ONasWUNISMjZsgMHDjBu3DiWLFnCzz//TGRkJP/+97/PlgcHB7Ns2TL69OmTr7+niIhIZmvWrKFDhw6kp6fzzTffcMUVV3i7SkXm28FUXuLjwR1ouFzOchHNmTPnbKDQp08f5syZA0CbNm24+OKLqVSpEn379mXZsmVUrlyZG264gY8//pjffvuN1NRUmjdvTvPmzVmyZAlPPvkk3333HVWrVj37/H/7298AaN68OVdddRVVqlShVq1aBAcHc+TIEay1PP300/z1r3/lpptuYteuXezduzfPei9dupQNGzawfv16Hn744bOtSQX1ww8/sHHjRq655hpatmzJjBkziM/0d+3du3ehnldEROTHH3+kY8eOBAcH891332W5e1KW+XbSzq+/zngcEJB1+ehROHwY0tKc5ZQUZ7ldO2e5Zs2s29epk+fhDh48yFdffcWGDRswxpCeno4xhltuueWciRXdy/fffz/PP/88TZo0oX///gA0atSI1atXs2jRIp566iluvvnmsy1P7g52fn5+WTrb+fn5kZaWRlxcHPv372f16tUEBAQQERHB6dOnsxw7OjqaTz75BHAi/MwuueQSateuzcaNG7noootwuVxny7I/jyfWWjp16nQ2iMyurCRQExER3/LNN9/QvXt3ateuzZdffkl4eLi3q1Rs8myZMsZMNcbsM8ZsyKE8yhiz7szP98aYFsVfTQ/GjnVaozJLTy9S36n58+dzzz33EB8fz/bt29m5cycNGzZk2bJl/PTTT2zbtg2Xy8XcuXNp3749AFdddRU7d+5k9uzZ9O3bF4Ddu3cTGhrKXXfdxeOPP372VmF+HD16lAsvvJCAgACWLl2apVXILSYmhjVr1pwTSAHs27ePbdu2ER4eTu3atdm3bx8HDx4kOTmZjz/++Ox2VapU4fjx4+fsf/XVV7N8+XL+/PNPAJKSkti8eXO+6y8iIpLd4sWL6dKlC2FhYXz77bflKpCC/LVMTQfeAP4vh/JtQAdr7WFjTFcgFriqeKqXixUrnNaozFJS4PvvC/2Uc+bMYcSIEVnW3XbbbUycOJG2bdsyYsQI1q9fz3XXXUfPnj3PbnPHHXewZs2asx2/169fz/Dhw/Hz8yMgIICJEyfmuw5RUVH06NGDyMhIWrZsSZMmTfK1X8eOHalUqRKpqam8+OKL1K5dG4BRo0Zx1VVX0bBhwyzP1a9fPwYNGkRISAgrVqw4u75WrVpMnz6dvn37kpycDMC4ceNo1KhRvs9BREQqtri4OKKjo9mxYwc1a9bk4MGD/PWvf+Xzzz+nVq1a3q5esTP56dxsjIkAPrbW5npz0xhTHdhgra2X13NGRkba7HmNNm3aRNOmTfOsj6/p3r07jz76KDfeeKO3q+J1ZfU1FBGR4pF9nj1wurJMnDiRgQMHerFmRWOMWW2tjfRUVtwd0O8DPs2lIgONMauMMav2799fzIcufUeOHKFRo0aEhIQokBIREcHzPHsul4vnn3/eSzUqecXWAd0Y0xEnmGqf0zbW2lic24BERkbm3STm46pVq6b+RCIiIpmU9Xn2CqNYWqaMMX8F3gZutdYeLMpz5ee2o/gmvXYiIhWbtZYqVap4LCsr8+wVRpGDKWNMA2ABcLe1tkjNNMHBwRw8eFAfymWQtZaDBw8SHBzs7aqIiIgXuFwuBg8ezLFjx86ZuLgszbNXGHne5jPGzAGuB2oaYxKA0UAAgLV2EjAKqAG8eSb3UlpOHbTyUr9+fRISEigP/akqouDgYOrXr+/taoiISClLTU2lX79+zJ49myeeeIK//vWvZ0fzNWjQgJiYmDI1PUxB5Ws0X0nwNJpPREREypZTp07Ru3dvPvroI1544YVzUgyVF7mN5vPtDOgiIiLis44dO8bf/vY3vv32W958800eeughb1fJKxRMiYiISIEdOHCArl278ssvvzBr1izuvPNOb1fJaxRMiYiISIHs2rWLm2++ma1bt/L+++/TvXt3b1fJqxRMiYiISL5t2bKFm266iQMHDvDpp59y/fXXe7tKXlfcGdBFRESknImLiyMiIgI/Pz8aNWrE/v37Wbp0qW8EUomJ0KED7NnjtSoomBIREZEcuefai4+Px1qLy+UiPT2d33//3dtVc4wdC8uWOb+9RKkRREREJEcRERHEx8efsz48PJzt27eXfoUyS0yEiy+G06chJAS2boU6dUrkUKU50bGIiIiUIz47196+fTB4MLhcznJ6utdapxRMiYiIiEcTJkzIcYo3r8+117s3LFwIKSnOckoKTJvmlb5TCqZEREQkC2sto0eP5uGHH6Z169aEhIRkKffaXHsrV8LJk87jCy+EgICs5V5qnVIwJSIiImelp6czePBgxowZQ//+/fnxxx+ZPHky4eHhGGMIDw8nNja29Ofa27oVrr4a/v1vZ3nzZkhNzbpNSgp8/33p1gt1QBcREZEzkpOTufvuu5k3bx5PPvkkL7zwAsYY71XI5YJ166BlS2d5zhzo1g3OP7/Uq6IO6CIiIpKr48eP061bN+bNm8err77Kiy++WPqBVPacUaNGQdu2kJDgLPft65VAKi/KgC4iIlLB7d+/n65du7JmzRpmzJjBPffc452KuHNGPf00TJ0KgwZBo0ZQr5536pNPus0nIiJSgW3fvp3OnTuzc+dO5s2bR7du3bxTkcw5oypVclqjSihnVGHoNp+IiIgAWaeGueiii2jVqhX79u3jiy++8F4gtW2b0yrlzhnl5+fVjOYFpWBKRESkgsg+NUxiYiJHjx5lxIgRXHPNNd6p1HvvwSWXwJQpGTmjUlO9ljOqMBRMiYiIVBDR0dEkJSVlWWetZeLEiaVbkdRU2L3bedy5M1xxxbnbeDGjeUGpA7qIiEgF4TNTw3TuDKdOwfLlcN55kJaW0Srl5qWcUYWhYEpERKQCsNZy/vnnc/To0XPKSmVqmIQEZ1SeMTBkCPj7O48Bfvml5I9fgnSbT0REpJxLS0vjwQcf5OjRo1SqVClLWYlNDZM5Z9Tq1XDppfDOO05Zz57Qo0dGMFXGKZgSEREpx5KSkujVqxeTJ0/m6aefZsaMGaUzNcyYMfDdd06/p5Yt4bHH4Npri/84PkB5pkRERMqpAwcO0L17d3766Sf+97//MXjw4NI5cGIihIU5nchDQpx59XwoZ1RhKM+UiIhIBbNt2zbatWvHmjVreO+990onkNq920m6OXaskysKytSovMJSMCUiIlLO/Pzzz7Rt25YDBw6wZMkSevbsWfIH3bkTLrvMub03bZqT/gCcUXllKGdUYSiYEhERKUc+//xzOnToQGBgIMuWLaN9+/YldzBrYdMm53FYGDz7LOzYkZHJ3K2ct07lGUwZY6YaY/YZYzbkUN7EGLPCGJNsjHm8+KsoIiIi+TFz5ky6devGxRdfzIoVK/jLX/5Ssgd86im46irYt89ZHj4cfv21TOeMKoz8tExNB7rkUn4IeAR4tTgqJCIiIvmTeZ696tWrc88993Dttdfy7bffUq9eveI9mDvVwapVsH+/s+6+++D116FmzYztfvnFabHK/lPGc0nlJs9gylr7LU7AlFP5PmvtSiC1OCsmIiIiOcs+z96RI0eoVKkS9957L1WrVi3+A44dC8uWQbt2MGqUs+6yy6Bfv4zO5hVUqZ69MWagMWaVMWbVfndUKyIiIgXmaZ699PR0Ro8eXbwHSkuDhQudTuTuvlD33Ve8xyjjSjWYstbGWmsjrbWRtWrVKs1Di4iIlCulNs/eyy9Dr15OJ3JwspZPm1a8xyjjKna7nIiISBn022+/4ZfDrbVimWdv1aqMUXq33goBARUq1UFBKZgSEREpQ7799lvatWtHaGgowcHBWcqKZZ69U6egSxcnzQHAhAnnzqFXzlMdFFR+UiPMAVYAjY0xCcaY+4wxg4wxg86U1zHGJAD/Ap45s835JVttERGRimf27Nl06tSJ2rVrs3btWt5+++3imWfv8GGYONEZdRcSAu+/D7GxTtmKFRUu1UFBaW4+ERERH2et5YUXXiA6OpoOHTqwcOFCqlevXvgnTEyEPn1g7lxnzrwJE2DIEFi7Fpo3L76KlyOam09ERKSMSk1N5YEHHiA6OpqoqCg+++yzogVS4Ez58t138OCDzvIDDyiQKgIFUyIiIj7q2LFjdO/enSlTpvDMM88wc+ZMgoKCivakiYkwfbpzS++TT5yO5IGBCqSKwN/bFRAREZFzJSQkcMstt7Bp0yamTJnCgAEDivaEa9bAa69BUFBGvqhKlZyO5BMmFLm+FZlapkRERHxE5ulhIiIi+OOPP1i0aFHRAymAnTvhww+dVil3h3KlOSgWCqZERER8QPbpYdLPJMnc555EuKCSk51M5e5Wp+7doWdP5/ZeZkpzUGQKpkRERHyAp+lhTp8+TXR0dO47uicgdrcupaU5vwMDnXWHDzvLxsDq1UpzUALUZ0pERMTL0tLSiI+P91iW5/Qw7gmIx46F666DJ590+kdVqwYff5w14eYvvxRbnSWDWqZERES86OjRo/To0SPH8lynh0lMzJiAeNo0qFkTrrgCTp50yrNnLpcSoWBKRETES7Zt28Y111zDkiVLuO+++wgNDc1Snuf0MKNHw+nTzuP0dFiwAN57D+rVK8FaS3YKpkRERLzg+++/56qrrmLXrl189tlnvP3228TGxuY9PUxamtPHKTERZs7MWK+ReV6jYEpERKSUxcXF0bFjR6pWrcqPP/7IDTfcAEBUVBTbt2/H5XKxfft2z/PsjR3rdDh/4omMfFFuGpnnFQqmRERESonL5WLUqFHcddddtG3blh9++IFGjRrltRO8+y5s2uQsP/igcytvwwaNzPMRCqZERERKwalTp+jbty9jx45lwIABfP7559SoUePcDbOnOjhyBO6/H956y1m+6CL429+ckXnWnvujEXulTsGUiIhICciczTwsLIxmzZoxb948XnnlFd5++20CAwM97zh2rDMJcZcuzvIFF8Dy5TB+fOlVXgpEeaZERESKmTubuTsJZ0JCAgDDhg3j8ccf97yTtRmpDqyFdeucW3tNm2oSYh+nlikREZFi5imbOcDChQs97xAfD23awMCBGZ3KAwLgjTdKsJZSXBRMiYiIFLOcspZnWe9ywe7dzuO6dZ3g6YsvNAlxGaRgSkREpBglJSUREhLisSxLNvPevaFzZyeoCgyEli3P3UGpDsoE9ZkSEREpJjt37uTvf/87SUlJBAQEUCM1lXeA3kBSSAgze/SA1FSnFap/f2cSYmudnVesUKqDMkrBlIiISDH4/vvv6dWrF0lJSXz00UccPXoU16BBtD9xgpfPO496gwdz7UsvQfv2TqvULbdkfQKlNCizFEyJiIgU0dSpUxk0aBDh4eEsXbqUpk2bOv2hkpMBuCc9HYYOhRtugJtu8nJtpbipz5SIiEghpaWlMXToUO677z6uv/56fvrpJ5q6M5qPG+fMowdO36dx4+Dmm8FPH73ljV5RERGRQjh48CBdunTh9ddf59FHH2XRokVUX7IELrkEfvstI18UaGReOadgSkREpIB+/fVX2rRpw3fffce7r7zCv0eMwN/fHy69FFq3hpgYTUJcgSiYEhERyUNcXBxX1q/PN8bQrGZNrrjiCpKSklj24YfcPnIkvPSSs2GrVrBggSYhrmDyDKaMMVONMfuMMRtyKDfGmNeNMX8aY9YZY1oXfzULLvM//pVhYcTFxXm7SiIiUga5p4bpv2sX7YHRBw8yMDmZp59+mis7d4bJk2HYsKw7aRLiCsVY9/3cnDYw5jrgBPB/1tpmHspvAYYAtwBXAa9Za6/K68CRkZF21apVhap0Xtz/+K8kJfEgMAl4IjSU2NhYoqKiSuSYIiJSPkVERJAcH89WIARIBQ4D14WF8VsOmc6l/DHGrLbWRnoqy7Nlylr7LXAol01uxQm0rLX2B6CaMaZu4apaPKKjozk/KYkBQCWgP1AlKYno6GhvVktERMqgavHxbCDjA9MFvA9sPjN5sUhx9JmqB+zMtJxwZt05jDEDjTGrjDGr9u/fXwyH9mzHjh2MBALPLPsBI8l5riQREZEsDhyAP/7gww8/JBmoBgSdKQoC7gKuqOfxo04qoOJI2mk8rPN479BaGwvEgnObrxiO7dEVF11E/127zkaKwTitU9P1jy8iInmxFnv11fwJ3LplC/9XpQqpx49TKdMmfsCsxo29VEHxNcXRMpUAhGVarg/sLobnLbS4Jk3OifD80T++iIhkkpgIHTo4uZ9WrYJHHwVrOXzkCKOrVePvW7bQv39/oho2JDjbrsFA44MHvVFr8UHF0TL1IfCwMeYdnA7oR621icXwvIXWyMM/eCBQb4PHAYkiIlIRPfccLFvm5H5q3RpmzmTTzTfT/eGH2blzJ69PnMiDDz6IMZ5uwIhkyDOYMsbMAa4HahpjEoDRQACAtXYSsAhnJN+fQBLOHTXvyjb0NOXECd5o1owxBw6wYtMmZ84kERGpuFauhNhYJ13BtGnw5JO8W6kS/W67jWrVqvHNN9/Qtm1bb9dSyog8UyOUlJJMjeBJQkICrVu3psEFF7B06VKq1PXqgEMRESltGzY4kw/ffDM89JATTLlc2MBAljdpwrXr1tG+fXvmzZtHnTp1vF1b8TG5pUaoMMEUwNLPPuOCLl04ERZGu/h4Nd2KiFQk118Pe/fCl1+SFhGBf2rq2aIkIGbAAEZPnEhgYGCOTyEVV5HyTJUnHTt3Zu/f/87InTv573//6+3qiIhIccncmdxt0SJo0QKOH3eW33oLli1j8z33kJYpkALnw/Ce+HgFUlIoFSqYAui0YAHVevZk+PDhfPfFF96ujoiIFIexY53O5MOGOa1PADVqQPXqGcuNG0ONGqR++63H0Xlp335bihWW8qRC3eZzO3r0KKObNuWxvXsJXL2a2i1beqUeIiJSDBIT4eKL4fRpZ/mRR+C11zxuevz4cc4//3yPZcYYXC5XSdVSyjjd5sumatWqPDxxIquNof+DD5KarblXRETKiFdega5dwR0E+fvDiRMeN92wYQNXXnlljk/VoEGDkqihVAAVMpgCuPTWWzk1cyaf/vQTTwwf7u3qiIhIfqSnw4oVGcubNsH69ZCS4iynpcGcOVn7TgH/93//R5s2bThy5AhPP/00oaGhWcpDQ0OJiYkp6dpLOVVhgymAvn37MnzQIDq89hrfP/KIt6sjIiLguTO524QJ0K6dE0QBBAU5rVGZpac7faiAU6dO8cADD3Dvvfdy5ZVX8ssvvxATE0NsbCzh4eEYYwgPDyc2NpaoqKgSPjEpryp0MAUw7tVXCataldlvvknL2rX5xhiuDAsjLi7O21UTEamY3J3Jx451Oo/36QNffumU9ekDc+fCJZc4yz/8kNEq5ZaSAt9/z5YtW2jXrh1vv/02I0aM4Msvv6TumRyDUVFRbN++HZfLxfbt2xVISZFU+GAqsHJlVowZw4T0dAbu20d7oF9CAgMHDlRAJSJS2hITnYzkLpfz+/RpZ1aLhASn/MIL4Y47wJ3C4JdfnCzm2X4WjBxJ69atiY+P56OPPuKFF17AP3sLlkgxqfDBFMCr//43dYD7gErAAKBKUhLR0dHerZiISEUzdiwkJzuP09Ph5Zfht9/g3ntz3CUuLo6IiAj8/PwIDw+na9eu3HbbbTRu3Jiff/6Z7t27l1LlpaJSMAXs2LGDkWT8MQwwEiA+3unMKCIiJeedd+Dqq2HnTqc1yp2yJyXFWXbnifIgLi6OgQMHEh8fj7WWHTt2sHjxYjp16sR3331HRERE6ZyDVGgKpoArLrqI/pyZvRkneVt/YIm/P/Ts6b2KiYiUB9k7lJ86BXFxcOiQs1y5MlSrBiNHZqQ4cMvUmdyT6OhokpKSzlm/efNmgoKCiukERHKnYAqIa9KE7LP0+QFpYWHw8MPOiuRkGDAA1q7NumFuo05ERCSjQ/mzzzrLmzfDXXfBggXOco8esHix8/6aQ2fynOzYsaNA60VKgoIpoNHBgx6nFkiOj2ete+bwjRvh/fczmpuPHYMjR7KOOhERkawydyifMsX54tmihZMrasCArNvm0JmcX37x+NRbt24lICDAY5kScEppUjAFHi/ghJ076Va3Lt26dWPXrl3QqhXs2gU33eTsM3EiXHRR1lEnap0SEXF89BHExDhfNDPfunN/8bz6avAr/EfQ3LlzadWqFZUqVTpncmIl4JTSpmAqB/Xr1+eTTz7h6NGjdO/enePHj0NISMbF37Wr8+3K/SaRkgL9+3uvwiIi3pSeDsuXZ3Qe/+YbmDzZ+aKZOTt5Eb94njx5kvvvv58+ffpw+eWXs3HjRqZOnaoEnOJVCqZy0aJFC+bNm8f69evp3bs3aZlH9tWqBWvWZLxJpKfDF19k7WApIlLeZO8n6g6epk2D9u1h3Tpn+dlns86Z55ZHh/LcrF+/niuvvJKpU6fy1FNP8c033xAREaEEnOJ1Cqby0KVLF958800+/fRThgwZgnW/cWRvugaoVMlZv2MH1KkDCxeWfoVFREqSu5/oE084rfPvv++s79nTyUzeqJGzfN55uWYnLwhrLRMnTuTKK6/k8OHDfPHFFzz//PM59pcSKW0KpvJh4MCBPPnkk0yaNIlXX33VWbliRe5vEv/4B1xxhfN440b45JOM4EsjAEWkrLEWXnkF3n7beS+bP99poQ8Jccpr1HAyk7uXocAdyj05dOgQt912G//85z/p2LEja9eu5cYbbyzmkxMpGgVT+fT8889zxx138MQTTzBv3rzc3yQaNHBGrbhHk7z5JvTtCydOOMsaASgiZcGhQ/Ddd85jY+DVV53bdOD8btwYunQp1kNmzmZep04dLrvsMj766CNeffVVPvnkEy688MJiPZ5IcTBnb1uVssjISLtq1SqvHLuwTp8+zY033sjq1av56quvaNeuXf52TE2F9euhdWunVap+feebXUgIbN3q3BIUESltiYkZEwe734dOnnSSaALcfbfTqr5nDxw8CBdf7MyV51bM72HubOaZk3AaY3juuecYOXJksRxDpLCMMauttZGeytQyVQDBwcF88MEHhIWFceutt7Jly5b87RgQ4ARSAM89l7E+PR3GjIHYWOeNSkSkNGVvJZ89G2rWhN27neURI5yBNQEBnvuJFqEzuSeesplba5kyZUqxHUOkJCiYKqCaNWuyaNEirLW0b9+esLAw/Pz8iIiIIC4uLvedExNhxoys6RSmToUHH4TPPnPWpaae+4YlIlLcEhOd95/MyTSvuAL++c+MbS6/3FlnTN79RIvIWkt8fLzHMmUzF1+nYKoQLrvsMgYPHsyePXtISEg4+yYwcODA3AMqT9/srHU6bbrnAJw502lKd38zBHVYF5HikZoKEyY4U7eMHZuR1sDdwtS4MYwf7yQkzq4YOpPn5NChQ/Tu3TvHcmUzF1+nYKqQZsyYcc66pKQkoqOjc94pp292mzdnjIAJD4cbb4S6dZ3lWbOgXz91WBeRvHn64rVpEyxd6jz293dG5MXFFXsyzcJasmQJzZs3Z+HChdxxxx2EhoZmKVc2cykL8hVMGWO6GGN+N8b8aYwZ4aG8ujFmoTFmnTHmJ2NMs+Kvqm8p1OSa+flmd+ONTpO7OTP18uuvw5IlGVPWfPNNxmgaEZHM3H2g/vWvjHWPPAIPPeQ8NgZWr4YqVUq8/1NeTp06xbBhw+jUqRNVq1blxx9/ZO7cucTGxiqbuZQ5eY7mM8ZUAjYDnYAEYCXQ11q7MdM2rwAnrLXPGWOaABOstbkmAimLo/kyi4iI8Hh/v0GDBjne9y+UQYMyvkEGBjpveI8/Di++6JRbmxF4iUjFlJQER49mHW23ZYuzvGEDVK8O9eplbN+qlTODQ3YtWxbLbbu8rFmzhqioKDZu3MiQIUN46aWXCMmcn0rEBxV1NF8b4E9r7VZrbQrwDnBrtm3+AnwJYK39DYgwxtQuQp19XkxMzDnN0QCNGzem2NJNuDusu5viU1KcLOtduzrLv/0Gl17quQOo+lmJlG/u95kPPnCCpccey2ht8veHl15yHjdrljWQghLt/5Rd5rxR4eHh9OnThzZt2nD48GEWL17M66+/rkBKyrz8BFP1gJ2ZlhPOrMtsLdALwBjTBggH6md/ImPMQGPMKmPMqv379xeuxj4iKioqS3N0gwYNuOWWW/jiiy8YPnx48QRUnjqsA7z7rvP71CknmIqIcJa//daZpf3ECSUGFSnrcvpC9McfcNllTv4ncFqT+vWDBQuy9oGaOdPrX6bceaPi4+Ox1rJjxw7mzp1Ly5YtWb9+PZ07d/Zq/USKS36CKU/3kLJHCi8C1Y0xa4AhwC9A2jk7WRtrrY201kbWqlWroHX1OZkn14yPj+fjjz/m4YcfZvz48YwaNaroB8hrKHKrVk5KBffIm6+/hv/8x8lZNW1a1iHPIlK2uL8QjR4NPXo4MymAM7NCs2ZOvydwBq1UqpTRUuVWyn2gPPGUNwpg79691KhRwws1EikZ/vnYJgEIy7RcH9ideQNr7TGgP4AxxgDbzvxUKMYYXnvtNZKTkxk3bhxBQUE888wzhX/Cgja5jxoFw4Y5ifYy57IaO9YZDr1tm/PG66dBnCI+behQmDzZuY5nzoR27TKCpaCgcydRL+EcUIWV04CcnTt3elwvUlbl51N1JXCZMaahMSYQ6AN8mHkDY0y1M2UA9wPfngmwKhw/Pz8mTZrE3XffzciRIzMmRi4tJ09mHfJsrbO8axdcdZWTINQt+5uv+lmJlIy8rq3Y2IwRdwCffprxhcg9B97gwTk/fyn2gcoPay2zZ8/G5DA4RnmjpLzJM5iy1qYBDwOfAZuAd621vxpjBhljBp3ZrCnwqzHmN6ArMLSkKlwW+Pn5MXXqVO644w6GDx/OG2+8UXoHz23Kh9deg/79nXX79zuzvM+enXVf9bMSKX7Zr61PPnGS9bpbm3budPJBuVxO4LVzZ9bWZS/lgCqMPXv20KtXL6KiomjYsCHBwcFZypU3Ssola61Xfq644gpb3qWkpNhbb73VAjY2NrZ0Dtqypafvp876zBISrB082Np165zljz+21s/P2TYkxNrExNKpr0h5t3u3tUFBzrUVHOxcW9OmWXv55dbu2eNs43JlbP/QQ9YGBma9fgMDrf3nP71S/fxyuVx29uzZ9oILLrBBQUH25ZdftmlpaXbWrFk2PDzcGmNseHi4nTVrlrerKlIowCqbQ0yjYKqEnT592nbt2tUaY+yMGTO8XZ2c9ehhrTEZb9ydOlnbq5e1R47kvM/u3dZed50CL6l48vrf/+MPa2+5xdoff3SCo4AA59oKCHCCoszBU3b5/ULkQxITE+3f//53C9irr77abtq0ydtVEil2uQVT6olcwoKCgnjvvfe44YYb6N+/P0OGDDmbcyVfkyOXhsREZ2Z49y2HlBRnZODWrXD++c66f/8bhgzJOmJItwWlonL/748Z4ywfOgTXX+9M0wJQrZoz4OOPP5xbdKmpzvrUVGd5796cn9vH+j/lxlrLnDlzuPzyy/n00095+eWXWbZsGU2aNPF21URKlYKpUhASEsIHH3zAZZddxhtvvHE250q+JkcuDZ76WRnjjCBydyBNTIQdOzKWBw/OGG1UhvpziBSatXD8uHMtuFOPTJ7s/O9Xrw4BAU6KAoCaNWHjRli+3OvTthSnzAk4w8LCuPLKK7nzzju57LLL+OWXXxg+fDiV3H8DkQpEwVQpqVy5ssd8K3lOjlwa8jOs+pVXnEzLbh9/nHW0Ubt2zgdMXjRiUHxNTv+TBw7AunUZy+3awd13Z/3yYa2zbIzTutunT9bn8NGUBYWRPQFnQkICq1evpk+fPixfvpymTZt6u4oiXqNgqhQlJCR4XJ/r5MiloaC3FRITYd++rKON4uOd9eBkZm/SBN55x1m2NmNb3RoUX+P+nxw6NOvo1v79swZH990HN9+cNfVIenruLbNl6JZdXnJKwLlixQq1RkmFp2CqFOWUWyUsLMzjep/l6bagv7+TywrgyBG4/HLnVgc4cwi60zC4b49MnarWKSl+BWn5/PJLuOuujP/JBQucAMrdv+npp2HSpIzt77/fmTS4HN22y6/09PQcJ3D3+pdBER+gYKoU5TQ5crNmzXB5moPPV+V166JuXXjvPbjpJmfZ3x9uv93JreM+z9OnnQ7t4HR0//TTjNnus9OtQcmvzC2f6enw+++QnOyULVgAYWFOjjVwOocvXJjxP+nnB337On2fANq2heuuy/r85ei2XX6tX7+ea665JsdyJeAUQakRSlvmnCsNGjSwnTp1soC99957bUpKirerV3J273Zy7GS+2eHOZ/Xyy87ygQPOtosWWTtypLWnTzvLDz3k5MDy8Tw7UgIKkn5j7dqMFAQhIdZOn+48/v57p/yHH6yNirI2Pt5ZTkjI+X9SbFJSkn3qqaesv7+/rVmzph00aJANDQ21OHOzWsCGhoYqb5RUGCjPlO9yuVz2ueees4Dt3r27PXnypLerVDJyS0R49Ki1y5dnbDtypLW1ajm5eHbvttbfP+sH3ebN1u7YkfOxlP+q/MgcSKelWbt+vbW7djll8fHWtm5t7fvvO8u9e2f93+rf30mOuXdvzs9dBpNjloYlS5bYSy65xAK2X79+9sCZLzpKwCkVmYKpMuDNN9+0xhjbvn17e+jQIW9Xp/gVNBFhcrLz+6GHrK1UKesHXffu1jZvnrHtf/9r7VtvZSwPGlTwliwFYCUvP3/jU6esPXzYebxzZ8ZrHxJi7e+/O49ffNEpP3nS2ptvtnbx4txbPnNSBpNjlrT9+/fbe+65xwL20ksvtV9++aW3qyTiMxRMlRFz5861AQEBtnnz5naX+9t3RZbTB+SiRdZ+/nnGdjfdZO3f/56xjzuTu/vDdPJka5cuzdjeU/Zp3UosmMIEn+6/8X33Zax74QVrM88MUL26M82Re/vsLUbz51u7davn51YrU4Fk73IwaNAgW6NGDevv72+jo6NtUlKSt6so4lMUTJUhX3zxha1cubKNiIiwf/zxh7er410F+YBMS8vYJ3tLVvXqWfepU8faZ57JWI6OzjhOSIi1GzfmPt2HW0EDCl9t/SpKYJT9tdi82drVqzOWX3rJ2piYrIGxn1/GsVq3trZfv4ztJ02y9quvCt7SpFamApk1a9Y5/Z8Ae8kll9j169d7u3oiPknBVBnz008/2Ro1atgLL7zQ/vzzz96ujvcU9AMypw/gLVusdd86TUuz9sknrf3gA2f5xAlnu8wBGDhBgLXWHjtmbefOTmuYtc6tpZkzrd2+PWtLS3p63udTmNavwgQ6Bd0nt3rFxzsdt90++MDa0aMz/s7+/tbecENGebdu1rZqlbHcp4+1//hH1sDY3z/jWDn93dTSVKLCw8PPCaQA26BBA29XTcRnKZgqgzZt2mTDwsJsUFCQrV27tjp85kdhPoA9BWABAdZ+9plTvmuXtZGR1i5Y4Cxv2uRs88YbWfeLjXXKf/3V2r/8xdqvv3aW//jDCba+/jpj+4AAa9esyTj+/PnWHjzoLO/fb+2yZU6QZ62zrzHWPvigs3zokDNqzd2nbO9eJ9hxLyckWPu3v2UER+vWOf3J3C13X3xh7bBhGS1v//ufs6078Bw2zNratTP+Pv/8p7U1amQsDxxobeXKGX/nSpWy9l9btcqZ3Devv7H6M3mVMcZjMGWM8XbVRHyWgqky6vXXXz/nTU9DkXNRmA/gggZgKSlOR+gBA7IGFFFRTvnvv1t7223WulsUv//e2osusrZnz6zH6dXLKf/wQ2d55Upnef58Z3ndOicIcQ/1Dwpygo9p05zlbduc7SdNcpbdfezGjs0asIwa5Tw+etQpf+kla6tUcTp6W+u0YLn7mAUGWtu1q7VPPJERbG3Y4Nx2c4uPL3hgpFYmn7Flyxbbo0cPj4EUYMPDw71dRRGfpWCqjMqpKV5veMWoMAFYQVtactv+2DEncHKnxNizx2kVO37cCULcwZQ7+Ni+3dr33stouYqPd25BuoOju+/OSCURGOi0bO3a5fl2WmFajAoTGKmVyetOnjxpR40aZYOCgmzlypVtnz59lDNKpIAUTJVRaor3UQUNKIrr9mNegU5B91FgVO65XC67cOHCs1/M+vbtaxMSEqy1yhklUlC5BVOaTsaH5TRNg7+/P9u3by/dykiGgk4pUpgpSDzNf5jXHHAF3acw9SpHE/eWN3FxcURERODn50dERATjx4+na9eu9OzZkypVqvD1118ze/Zs6tWrB0BUVBTbt2/H5XKxfft2oqKivHwGImWXcYKt0hcZGWlXrVrllWOXFXFxcQwcODDLTO1BQUH4+fkRGhrKu+++yw033ODFGkqJadUK1qw5d33LljkHLoXZR8oFT+8VAMHBwbz44ov885//JMA956CIFIoxZrW1NtJTmVqmfFhUVBSxsbGEh4djjCE8PJwpU6awdu1aLrzwQm6++WZee+01vBUQSwkqTAuQWo0qrOjo6HMCKYCaNWsydOhQBVIiJUwtU2XUsWPHuOeee/jggw+49957mTRpEsHBwd6uloh4gZ+fn8cvVcYYXNlv/YpIoahlqhw6//zzWbBgAc8++ywzZszguuuuIyEhwdvVEpFStG3bNvr06ZNj63RO/S5FpHgpmCrD/Pz8GD16NO+//z6bNm0iMjKS5cuXe7taIlLCDh8+zPDhw2nSpAkffvghPXv2JCQkJMs2oaGhxMTEeKmGIhWLgqly4NZbb+XHH3/k/PPPp2PHjgwYMCDLqJ64uDhvV1FEikFKSgqvv/46l156KePHjycqKoo//viDBQsWMHny5Cz9K2NjYzVCT6SUqM9UOXLkyBE6dOjAunXrsqwPDQ3VG6tIGWat5f333+eJJ57gzz//5KabbuLVV1+lRYsW3q6aSIVR5D5TxpguxpjfjTF/GmNGeCivaoz5yBiz1hjzqzGmf1ErLQVXrVo1jhw5cs76pKQkoqOjS79CIlIomXNG1a1bl6ZNm9KrVy8CAwNZtGgRn3/+uQIpER/in9cGxphKwASgE5AArDTGfGit3Zhps8HARmttD2NMLeB3Y0yctTbFw1NKCdq5c6fH9Tt27CjlmohIYWTPGbVnzx727NnDgAEDeOutt/D3z/NtW0RKWX5aptoAf1prt54Jjt4Bbs22jQWqGGMMcB5wCEgr1ppKvuQ0eic4OJh9+/aVcm1EpKCefPJJjzmjvvzySwVSIj4qP8FUPSBzc0fCmXWZvQE0BXYD64Gh1tpzkpsYYwYaY1YZY1bt37+/kFWW3MTExBAaGpplXUBAACkpKTRv3pyPP/7YSzUTkdzs37+fxx9/nF27dnksV+uyiO/KTzBlPKzL3mu9M7AGuAhoCbxhjDn/nJ2sjbXWRlprI2vVqlXAqkp+eMqaPm3aNNasWUPdunXp0aMHgwYN4uTJk96uqojgpDl45plnaNiwIf/5z3+oXLmyx+2UM0rEd+UnmEoAwjIt18dpgcqsP7DgzMTKfwLbgCbFU0UpKE8TmDZr1owff/yR4cOHExsbS6tWrfjpp5+8XVWRCuv48ePExMRw8cUXExMTQ/fu3fn111956623zmldVs4oEd+Wn2BqJXCZMaahMSYQ6AN8mG2bHcCNAMaY2kBjYGtxVlSKLigoiJdffpmlS5eSnJxMu3btGDNmDGlp6t4mUlIyj8yLiIhg2rRpjB8/nosvvphnnnmGDh06sHbtWt555x2aNGnisXVZqU1EfFu+8kwZY24B/gtUAqZaa2OMMYMArLWTjDEXAdOBuji3BV+01s7K7TmVZ8q7jhw5wsMPP0xcXBxXX301d9xxB6+99ho7duygQYMGxMTE6M1bpIiyj8zLrHPnzowZM4Y2bdp4oWYiUlC55ZlS0s4K7p133mHAgAGcOnUqy3ol+hQpuoiICOLj489ZX7t2bfbs2eOFGolIYSmYklzVr1/f4wii8PBwtm/fXvoVEikHjh07RtWqVT2WGWNwuc4Z8CwiPqzIGdClfNu9O/t4AoeGYosU3P79+3nmmWcIDw/PcRuNzBMpXxRMSY5v7H5+fsyfPx9vtV6KlCU7duxg6NChhIeH8/zzz3PDDTcwZswYjcwTqQAUTInHRJ9BQUHUq1eP22+/nc6dO/P77797qXYiviX76LxXXnmF/v37c8kll/Dmm2/Su3dvNm7cyHvvvcfIkSM1Mk+kAlCfKQGcD4jo6Ogso/l69+7NpEmTeOaZZ0hKSuKxxx7jmWeeyTGpoEh5l9PovICAAB566CEee+wx3cITKafUAV2KZO/evYwYMYLp06cTFhbGf/7zH06fPn1O8KVv21LehYeHe+xLWL9+/RwnGReR8kHBlBSL5cuXM3jwYNauXYufn1+W0UhKpSDl2dGjR5k+fTrDhg3zWK7ReSLln4IpKTZpaWlceOGFHD58+JwypVKQ8ub333/njTfeYPr06Zw4cYKgoCCSk5PP2U7/+yLln1IjSLHx9/fnyJEjHsuUSkHKA5fLxaJFi+jSpQtNmjQhNjaWXr16sXLlSqZMmaLReSJyDgVTUmC5dbAdN24cx44dK8XaiBRe5pF5DRo04O6776Zx48Z069aNdevWMXbsWHbu3MmMGTOIjIzUvHki4pFu80mBeRrRFBwcTNOmTfnll1+44IILePzxxxkyZAjnnXeeF2sqkrOcRuZddtlljBkzhl69ehEYGOil2omIr9FtPilWnr6dv/322/z888/89NNPXH311Tz99NM0bNiQl19+mZMnT3q7yiJZHDt2jEceecTjBMQpKSn06dNHgZSI5JtapqRE/Pjjj4wePZrPPvuMWrVq8eSTT3LBBRfw3HPPKZ2CeIW1lm+//ZapU6cyf/58j4EUaGSeiHimlikpdVdddRWLFy9m+fLltGjRgscff5wBAwYQHx+PtZb4+HgGDhxIXFyct6sq5UT2zOTu/62dO3cybtw4Lr30Uq6//nref/997rrrLurUqePxeZR0U0QKSi1TUirq1KnD3r17z1mvIeVSHDz1fwoMDKRRo0b8+uuvWGvp2LEjAwYMoFevXoSGhnrcR/nSRCQnapkSr9u3b5/H9fHx8Xz66ae6rSJFEh0dfc5tu5SUFDZt2sTIkSPZunUrX331FXfdddfZ1AYamScixUUtU1IqIiIiiI+PP2e9O5P6xRdfzKBBgxgwYAA1atTwQg2lLHK5XKxYsYL27dt7LFf/JxEpLmqZEq+LiYnxmOxw6tSpzJkzh3r16vHEE09Qr149+vXrx8qVK4Gc+8FIxeVyuVi2bBlDhw6lQYMGOQZSoP5PIlI61DIlpSYuLi7XyZHXr1/PxIkTmTlzJidOnKBhw4bs2rWLlJSUs9uoT0vFkfn/JSwsjHvvvZcjR47w3nvvsXv3boKCgujSpQu33347p06dYujQoer/JCIlRnPzSZly7NgxZs6cybBhw0hLSzunXJ3Wy7+cEmr6+/vTo0cPbr/9drp168b555+fZZ/cgnURkaJQMCVlkp+fHzn9fy5evJgbbriBgICAUq6VlKT4+Hg+++wzhg0bxqlTp84pDwsL0xyQIuIV6jMlZVJO/V2MMXTp0oXatWvTv39/Fi1adPZWoPpY+aacXpdTp06xePFiHn30UZo2bUpERAQPPvigx0AKICEhoTSrLSKSL2qZEp+VUx6gCRMmUKNGDebNm8eHH37I0aNHqVq1Ks2aNWPVqlUkJydn2V79ZrzL0+sYEBBA48aN+fPPPzl9+jTBwcF06NCBzp0706VLF7p06eKxBUq3eEXEW3SbT8qsvPrBJCcn8+WXXzJ//nxmzJjhcRh8gwYNPKZlkJLncrmoX78+iYmJ55T5+/szZMgQOnfuzHXXXUdISMjZMiXUFBFfk1swhbXWKz9XXHGFFSlOxhgLePzp27evnTRpkt20aZN1uVxZ9ps1a5YNDw+3xhgbHh5uZ82a5aUzKBty+3sdPXrUfv755/a5556znTt3ttWqVcvxNTHGFPo4IiKlDVhlc4hp8tUyZYzpArwGVALetta+mK18OOD+uugPNAVqWWsP5fScapmS4pZTYtDQ0FCqVq16tnWkdu3adOjQgQ4dOnD8+HHGjBmjFpB8ymnalmuuuYYDBw6wYcMG543FGJo1a0bbtm157733OHjw4DnPpVt2IlKWFKllCieA2gJcDAQCa4G/5LJ9D+CrvJ5XLVNS3GbNmmVDQ0OztH6EhobaWbNmWZfLZTdv3mxjY2NtVFSUrVevXo4tJoBt0KBBnscqD60m+T2PQ4cO2WXLltkLLrggx1amzp072+eee85+/vnn9siRI1mOkdPrIiJSVpBLy1R+gqm2wGeZlp8Cnspl+9nAA3k9r4IpKQn5DQ5cLpf9888/cw2o2rVrZwcOHGhfe+01u2TJErtnzx7rcrl8OjgoSJDn6TxCQkLsqFGj7FtvvWUfeeQRe9NNN9m6devm+ndyB1PFVS8REV9U1GDqHzi39tzLdwNv5LBtKHAIuCCv51UwJb4gPDzcY3Bw3nnn2WuvvfaclpgaNWrYoKCgArdmFSaYKOg+eQV56enpds+ePXb16tX2ww8/zLGVKfO+kZGR9t5777UvvfSS/fjjj3Ns0QsPDy/Q311EpKzJLZjy93jvLyvjYV1OHa16AMttDn2ljDEDgYGgObPEN8TExHgcNTZp0iSioqKw1rJnzx5+/fXXsz+TJ0/2+Fw7duygdu3a1KlTh7p161KnTh3q1KnDrl27ePfdd8/mwoqPj+eBBx7gxIkT3HnnnQQGBhIYGIgxGZda9r5J8fHxDBw4EIDbb7+dkydPnv1JSkri5MmT/Otf/zonY3hSUhL33Xcf0dHR7N69m9TU1Dz/JsYYtm7dSoMGDfDzy5qK7qWXXvL494qJicnzeUVEyqs8O6AbY9oCz1prO59ZfgrAWvuCh20XAvOstbPzOrA6oIuvKOg0JDl1dK9atSq9e/cmMTGRPXv2nP3taUocTwICAggMDCQoKIgjR454TPNQWHfddRf16tWjXr161K9fn3r16tGrVy927dp1zrZ5dQzXtC0iUhEVKc+UMcYf2AzcCOwCVgJ3Wmt/zbZdVWAbEGatPZlXpRRMSVlVkBxILpcLf3//HKfFeeWVV0hJSSE5OZmUlJSzjydMmJDj8ceNG0doaCiVK1fO8nP33Xezd+/ec7bPKThSLicRkfwrcp4p4BacgGoLEH1m3SBgUKZt+gHv5Of5rPpMSRlXkP5MOfXLyq2fUWH2KUzHeHUMFxHJH4qaZ6okqGVKKorCtAAVttVIt+BEREqGppMR8bLCBDkKjEREfIeCKREREZEiyC2Y8vO0UkRERETyR8GUiIiISBEomBIREREpAgVTIiIiIkWgYEpERESkCLw2ms8Ysx84d06O4lcTOFAKx/FFOveKqyKff0U+d6jY569zr7hK4/zDrbW1PBV4LZgqLcaYVTkNZSzvdO4V89yhYp9/RT53qNjnr3OvmOcO3j9/3eYTERERKQIFUyIiIiJFUBGCqVhvV8CLdO4VV0U+/4p87lCxz1/nXnF59fzLfZ8pERERkZJUEVqmREREREqMgikRERGRIiizwZQxposx5ndjzJ/GmBEeyo0x5vUz5euMMa3zu6+vy8e5R50553XGmO+NMS0ylW03xqw3xqwxxqwq3ZoXj3yc//XGmKNnznGNMWZUfvf1dfk49+GZznuDMSbdGHPBmbIy/dobY6YaY/YZYzbkUF5ur3nI1/mX2+s+H+denq/5vM69PF/zYcaYpcaYTcaYX40xQz1s4xvXvbW2zP0AlYAtwMVAILAW+Eu2bW4BPgUMcDXwY3739eWffJ57O6D6mcdd3ed+Znk7UNPb51HC53898HFh9vXln4LWH+gBfFWOXvvrgNbAhhzKy+U1X4DzL8/XfV7nXi6v+fyce7Zty9s1XxdofeZxFWCzr37Wl9WWqTbAn9bardbaFOAd4NZs29wK/J91/ABUM8bUzee+vizP+ltrv7fWHj6z+ANQv5TrWJKK8vqV+9c+m77AnFKpWSmw1n4LHMplk/J6zQN5n395vu7z8drnpMy/9gU89/J2zSdaa38+8/g4sAmol20zn7juy2owVQ/YmWk5gXP/wDltk599fVlB638fTtTuZoHPjTGrjTEDS6B+JS2/59/WGLPWGPOpMebyAu7rq/Jdf2NMKNAFeC/T6rL+2uelvF7zhVHervv8KI/XfL6V92veGBMBtAJ+zFbkE9e9f0k9cQkzHtZlz/GQ0zb52deX5bv+xpiOOG+q7TOtvsZau9sYcyHwhTHmtzPffMqK/Jz/zzhzKJ0wxtwCvA9cls99fVlB6t8DWG6tzfyNtqy/9nkpr9d8gZTT6z4v5fWaL4hye80bY87DCRKHWWuPZS/2sEupX/dltWUqAQjLtFwf2J3PbfKzry/LV/2NMX8F3gZutdYedK+31u4+83sfsBCnKbQsyfP8rbXHrLUnzjxeBAQYY2rmZ18fV5D69yFbc385eO3zUl6v+Xwrx9d9rsrxNV8Q5fKaN8YE4ARScdbaBR428Y3rviQ7j5XUD06L2lagIRkdyy7Ptk03snZK+ym/+/ryTz7PvQHwJ9Au2/rKQJVMj78Hunj7nErg/OuQkZC2DbDjzP9BuX/tz2xXFaePReXy9NqfqXsEOXdCLpfXfAHOv9xe9/k493J5zefn3M+Ul8tr/sxr+H/Af3PZxieu+zJ5m89am2aMeRj4DKfH/lRr7a/GmEFnyicBi3B6+f8JJAH9c9vXC6dRKPk891FADeBNYwxAmnVm064NLDyzzh+Yba1d7IXTKLR8nv8/gIeMMWnAKaCPda6uivDaA/QEPrfWnsy0e5l/7Y0xc3BGbdU0xiQAo4EAKN/XvFs+zr/cXvf5OPdyec1Dvs4dyuk1D1wD3A2sN8asObPuaZwvDj513Ws6GREREZEiKKt9pkRERER8goIpERERkSJQMCUiIiJSBAqmRERERIpAwZSIiIhIESiYEhGfZ4ypZoz555nHFxlj5nu7TiIibkqNICI+78y8XB9ba5t5uy4iItmVyaSdIlLhvAhcciZx3x9AU2ttM2NMP+DvOEn5mgHjcbId3w0kA7dYaw8ZYy4BJgC1cBL7PWCt/a20T0JEyifd5hORsmAEsMVa2xIYnq2sGXAnzjQiMUCStbYVsAK458w2scAQa+0VwOPAm6VRaRGpGNQyJSJl3VJr7XHguDHmKPDRmfXrgb+emXG+HTDvzNQaAEGlX00RKa8UTIlIWZec6bEr07IL5z3ODzhyplVLRKTY6TafiJQFx4EqhdnRWnsM2GaMuR3AOFoUZ+VEpGJTMCUiPs9aexBYbozZALxSiKeIAu4zxqwFfgVuLc76iUjFptQIIiIiIkWglikRERGRIlAwJSIiIlIECqZEREREikDBlIiIiEgRKJgSERERKQIFUyIiIiJFoGBKREREpAj+H5ZMggdRMRVgAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "tags": [], "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "pK6dv1hQjODA" }, "source": [ "## Consistency \n", "To prove that the Abysmal-Butler method does not satisfy the consistency condition,\n", "\\begin{equation} \\tau_{i}(h)=\\frac{y_{i+1}-y_{i}}{h}-\\frac{1}{2}[4(f(t_i,y_i)-3f(t_{i-1},y_{i-1})]. \\end{equation} \n", "As $h \\rightarrow 0$ \n", "\\begin{equation} \\frac{1}{2}[4f(t_i,y_i)-3f(t_{i-1},y_{i-1})] \\rightarrow \\frac{f(t_i,y_i)}{2}.\\end{equation} \n", "While as $h \\rightarrow 0$ \n", "\\begin{equation} \\frac{y_{i+1}-y_{i}}{h} \\rightarrow y^{'}=f(t_i,y_i).\\end{equation} \n", "Hence as $h \\rightarrow 0$ \\begin{equation} \\frac{y_{i+1}-y_{i}}{h}-\\frac{1}{2}[4(f(t_i,y_i)-3f(t_{i-1},y_{i-1})]\\rightarrow f(t_i,y_i)-\\frac{f(t_i,y_i)}{2}=\\frac{f(t_i,y_i)}{2},\\end{equation} \n", "which violates the consistency condition (inconsistent).\n" ] }, { "cell_type": "code", "metadata": { "id": "5WrT05UxjODB", "outputId": "4d87cab9-8cc2-4bd1-964e-50070dc1be6b" }, "source": [ "d = {'time': t[0:5], 'Abysmal Butler': w[0:5],'Exact':y[0:5],'Error':np.abs(y[0:5]-w[0:5])}\n", "df = pd.DataFrame(data=d)\n", "df" ], "execution_count": null, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
timeAbysmal ButlerExactError
00.001.0000001.0000000.000000
10.050.9524590.9524590.000000
20.100.9372130.9096750.027538
30.150.9211760.8714160.049760
40.200.9068490.8374620.069388
\n", "
" ], "text/plain": [ " time Abysmal Butler Exact Error\n", "0 0.00 1.000000 1.000000 0.000000\n", "1 0.05 0.952459 0.952459 0.000000\n", "2 0.10 0.937213 0.909675 0.027538\n", "3 0.15 0.921176 0.871416 0.049760\n", "4 0.20 0.906849 0.837462 0.069388" ] }, "metadata": { "tags": [] }, "execution_count": 11 } ] }, { "cell_type": "code", "metadata": { "id": "7pdVDKuPjODD" }, "source": [ "" ], "execution_count": null, "outputs": [] } ] }